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Preface 


This  report  is  the  culmination  of  my  efforts  to 
understand  selected  physical  phenomena  inherent  in  a 
thermonuclear  plasma.  An  incomplete  computer  simulation 
of  inertially  confined  fusion,  program  MOXNEX,  is  the 
central  vehicle  in  the  study.  It  is  hoped  that  the  modest 
development  and  documentation  contained  here  may  contribute 
to  the  understanding  of  others  participating  in  programs 
in  Nuclear  Engineering  and  Physics  at  the  Air  Force 
Institute  of  Technology. 

I  would  like  to  thank  Doctor  George  Nickel,  not  only 
for  originating  the  project,  but  for  extending  it  for 
further  study.  I  would  like  to  thank  Major  Michael  Stamm 
and  especially  my  advisor,  Lieutenant  Colonel  William  Bailey, 
for  patient  consultation,  advice,  guidance  and  leadership. 

I  would  also  like  to  thank  Ms.  Sharon  Gabriel  for  manuscript 
preparation.  I  thank  Miss  Patricia  Horton  for  patience, 
understanding  and  support.  Finally,  and  most  importantly, 

I  would  like  to  thank  my  parents,  Joseph  and  Florence 
DeBruyne,  for  continued  faith  and  encouragement  in  this  and 
other  endeavors. 
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Abstract 


Energy  transfer  processes  in  a  thermonuclear  plasma 
including  Coulombic  phenomena,  bremsstrahlung ,  and  Compton 
scatter  are  critically  reviewed.  These  processes  are 
incorporated  in  a  three  temperature,  inertially  confined 
fusion  computer  simulation  which  uses  a  one  dimensional, 
spherical,  Lagrangian  hydrodynamics  scheme.  The  computer 
code,  still  in  the  validation  phase,  uses  separate  subroutines  to 
model  hydrodynamics ,  thermonuclear  bum,  neutron  heating,  alpha  par¬ 
ticle  heating,  and  energy  transfer  processes  in  CGS  units. 
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one  dimensional  analysis 

OF  INERTIALLY  CONFINED  PLASMAS 

I.  Introduction 

Background 

The  scarcity  of  energy  resources  is  a  growing  national 
concern.  Nuclear  energy  has  helped  to  ensure  an  independent 
and  balanced  energy  supply  system  since  the  early  1960’s 
through  the  use  of  fission  reactors.  An  energy  source  yet 
to  contribute  to  civil  energy  reserves  is  nuclear  fusion. 

Nuclear  fission  is  the  exothermic  breakup  of  certain 
heavy  elements  into  lighter  components  and  is  the  driving 
mechanism  for  the  nuclear  reactors  in  service  today.  Nuclear 
fusion  is  the  exothermic  combination  of  certain  light  nuclei 
into  a  heavier  nucleus.  Two  concepts  show  promise  in 
exploiting  nuclear  fusion  as  an  energy  resource.  Magnetic 
fusion  uses  a  magnetic  field  to  confine  a  thermonuclear 
plasma.  Inertial  confinement  fusion  bombards  a  fuel  pellet 
with  a  concentrated  energy  pulse  creating  a  thermonuclear 
plasma  by  the  resulting  compression  and  heating.  Both  con¬ 
cepts  warrant  further  study. 

Three  technologies  are  being  researched  to  supply  the 
concentrated  energy  pulse  for  an  inertial  confinement 
scheme,  namely:  electron  beams,  ion  beams,  and  laser  beams. 
These  energy  sources,  known  as  drivers,  couple  energy  into 
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a  fuel  target  only  micrometers  in  diameter.  The  surface 
layer  of  the  target  absorbs  much  more  energy  than  the 
target  interior.  As  a  result,  the  outer  layer  vaporizes 
and  expands  outward  extremely  rapidly.  This  outward  pres¬ 
sure  must  be  balanced  by  inward  pressure  which  is  manifest 
in  a  shock  wave.  The  shock  wave  compresses  the  fuel  and 
the  compression  action  raises  the  temperature.  The  result¬ 
ing  thermonuclear  burn  will  free  large  amounts  of  energy. 

In  fact,  this  is  a  nuclear  explosion;  it  is  a  miniature 
nuclear  bomb.  Very  careful  design  of  the  target  using 
sophisticated  physics  optimizes  the  pellet  compression  and 
attains  a  maximum  yield  from  the  fuel  mixture  (see  figure  1). 

The  U.S.  government  has  identified  nuclear  fusion  as 
an  inexhaustible  energy  source  and  sees  its  commercial  use 
by  2020  (Ref  19:133-138).  In  view  of  the  continuing  need 
for  energy,  the  government  has  budgeted  $593.7  million 
dollars  for  fiscal  year  1981  (including  Reagan  cuts)  and 
$650.3  million  dollars  for  fiscal  year  1982  for  fusion  energy 
research.  Of  these  totals,  inertial  confinement  fusion  will 
account  for  $199.6  million  and  $190.2  million  for  operating, 
capital  equipment,  and  construction  for  fiscal  years  1981 
and  1982,  respectively  (Ref  20:55-57). 

During  the  Winter  Quarter  1980  at  the  Air  Force  Institute 
of  Technology,  Dr.  George  H.  Nickel,  assisted  by  class  members, 
developed  the  computer  code  NOXNEX  to  model  inertially  confined 
fusion  in  DT  microspheres  while  instructing  a  graduate  class  in 
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nuclear  fusion.  The  project  was  founded  primarily  on  two 
prior  studies,  one  completed  by  Keith  A.  Brueckner  and 
Siebe  Jorna  at  KMS  Fusion  in  1974  (Ref  4) ,  and  the  other 
done  by  a  team  headed  by  Gary  S.  Fraley  at  Los  Alamos 
Scientific  Laboratory,  also  in  1974  (Ref  9).  Prior  to  his 
retirement  from  the  Air  Force  and  subsequent  assumption  of 
a  position  at  the  Los  Alamos  National  Laboratory,  Dr.  Nickel 
entrusted  the  code  to  Dr.  William  F.  Bailey  for  further 
study  and  development. 

Goals  and  Discussion 

The  goals  of  this  project  are  to  develop  a  computer 
code,  MOXNEX,  to  model  an  inertially  confined  fusion  plasma, 
validate  the  code,  and  prepare  the  code  for  further  studies 
in  Nuclear  Engineering  and  Physics  at  the  Air  Force  Institute 
of  Technology.  Though  the  code  has  yet  to  model  fusion, 
progress  has  been  made  toward  all  three  goals. 

As  previously  mentioned.  Reference  9,  published  in 
The  Physics  of  Fluids,  was  a  major  source  for  initial  code 
development.  Upon  critical  examination,  however,  inconsis¬ 
tencies  became  apparent  in  the  equations  of  Reference  9 
in  the  area  of  radiation-electron  energy  transfer.  These  are 
mentioned  in  the  text  and  appendices.  Dr.  Nickel's  efforts  to 
resolve  these  inconsistencies  through  conrnunication  with  the  authors 
were  unfruitful.  Development  of  consistent  formulation,  then, 
was  necessary.  The  development  of  the  radiation-electron 
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energy  transfer  equations  is  one  of  the  major  topics  of 
this  report. 

Additionally,  simplifications  were  included  in  the 
coding  to  expedite  code  construction.  One  such  simplifi¬ 
cation  was  the  bremsstrahlung  Gaunt  factor  being  set  equal 
to  1.  A  complicated  integral  was  evaluated  numerically 
and  the  resulting  values  were  fit  with  an  exponential 
function. 

Validation  is  required  to  enable  its  use  as  a  vehicle 
for  further  study.  Confidence  in  the  MOXNEX  code  will  allow 
its  application  to  numerous  areas  of  education  and  research. 
Code  hydrodynamics  were  examined  to  provide  some  portion  of 
this  confidence. 

Finally,  code  documentation  is  necessary  to  provide 
user  utility.  Under  Dr.  Nickel's  tutelage,  each  member  of 
the  Fusion  Engineering  and  Devices  class,  taught  during  the 
Winter  Quarter  of  1980,  contributed  to  the  construction  of 
the  MOXNEX  code.  Because  of  the  diversity  of  the  originators, 
the  code  is  somewhat  disjointed  and  is  not  strongly  supported. 
Notes  from  each  of  the  class  members  and  from  Dr.  Nickel 
are  assimilated  and  presented  in  this  report. 

The  presence  of  an  inertial  fusion  code  at  the  Air 
Force  Institute  of  Technology  provides  Physics  Department 
programs  with  a  valuable  tool.  Classes  in  nuclear  engineer¬ 
ing,  engineering  physics  and  electro-optics  may  easily  find 
uses  for  the  MOXNEX  code  and  its  results.  Classes  in 
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nuclear  explosives,  fusion  engineering,  plasma  physics 
laser  studies,  and  others  may  be  directly  supported  by 
code  use. 

Additionally,  further  code  development  might  also  be 
undertaken  as  design  study  or  thesis  projects.  Extension  of 
the  code  to  include  driver  input  would  make  study  of  laser - 
plasma  coupling  possible. 

With  this  in  mind,  much  work  is  required  for  completion 
of  the  MOXNEX  code.  It  is  hoped  that  it  can  be  a  valuable 
vehicle  for  the  Department  of  Physics  at  the  Air  Force 
Institute  of  Technology. 

Scope 

This  analysis  is  concerned  with  the  thermonuclear 
burn  and  related  hydrodynamic  and  energy  transfer  processes, 
code  validation,  and  documentation  of  the  coding.  Mechanisms 
for  the  initial  delivery  of  energy  were  not  to  be  considered, 
nor  was  the  coupling  of  energy  to  the  plasma  investigated. 

Code  Synopsis  with  Assumptions 

The  three  temperature  MOXNEX  code  models  inertially 
confined  fusion  in  deuterium-tritium  microspheres  using 
separate  ion,  electron,  and  radiation  temperatures.  The  ions 
and  electrons  are  modelled  by  a  Maxwellian  energy  distribution 
and  the  photons  are  modelled  by  a  Planckian  distribution. 

The  program  is  written  in  Fortran  V  using  the  central  program 
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to  call  each  of  seven  primary  subroutines.  A  general 
flowchart  is  seen  in  Figure  2.  The  subroutines,  tasks 
performed,  and  assumptions  are  discussed  briefly  in  the 
following  paragraphs.  Constant  temperatures  and  number 
densities  are  assumed  at  the  time  of  each  subroutine  call. 
Appendix  A  is  a  program  listing  and  Appendix  B  is  a  glossary 
of  variables. 

1.  Subroutine  GDATA.  This  subroutine  sets  initial 
conditions  by  initializing  variables  and  is  called  one 
time  prior  to  execution  of  physical  processes.  Energy  is 
deposited  in  the  microsphere  by  initializing  ion,  electron, 
and  radiation  temperatures  to  selected  levels  at  selected 
cell  locations.  Compression  is  also  set  as  an  initial 
condition  by  specifying  an  initial  density. 

2.  Subroutine  HYDRO.  Subroutine  HYDRO  is  a  one¬ 
dimensional  spherical  Lagrangian  hydrodynamics  code.  It 
accounts  for  the  hydrodynamics  of  the  cells  and  the  change 
in  temperatures  due  to  PdV  work. 

Assumptions  in  HYDRO  include  (1)  all  materials  in 
the  microsphere  are  an  ideal  gas,  (2)  pseudo-viscous  pressure  is 
due  only  to  the  ions,  (3)  electron  degeneracy  occurs  only 
in  the  deuterium-tritium  fuel  regions  and  can  be  approxi¬ 
mated  with  an  effective  temperature  to  account  for  additional 
pressure,  (4)  total  system  PdV  work  changes  can  be  divided 
between  the  electron,  ion,  and  radiation  species  based  on 
the  pressure  of  each  component.  Viscous  pressure  is  set 
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Figure  2.  General  flowchart  for  Program  MOXNEX 


to  zero  when  the  cell  has  expanded  during  the  previous  time 
step . 

3.  Subroutine  TBURN.  The  TBURN  subroutine  provides 
reaction  rates  of  the  deuterium-tr itium  fusion  reaction. 
Deuterium-deuterium  reactions  are  not  considered.  Tritium 
created  by  deuterium-deuterium  reactions  is  not  included. 

If  the  ion  temperature  is  below  a  specified  ignition  temper¬ 
ature  of  1  kev,  a  figure  based  on  bremsstrahlung  power  lost, 
no  reactions  are  computed  in  any  cell.  The  number  of  deu¬ 
terium  and  tritium  ions  burned  are  subtracted  from  the  totals 
in  each  cell. 

The  hydrogen  fuel  mix  is  equal  amounts  of  deuterium  and 
tritium  at  all  times  in  TBURN.  Therefore,  the  hydrogen  number 
density  is  used  in  modelling  the  reaction  rate  density. 

The  Maxwellian  velocity  distribution  weighted  average  of 
cross  section  for  the  deuterium- tritium  reaction  is  modelled 
in  two  ranges  of  ion  temperature.  Curve  fits  to  data  above  10 
kev  below  10  kev  ion  temperature  are  used  for  the  model. 

4.  Subroutine  NUHEAT.  Subroutine  NUHEAT  calculates 
total  neutron  heating  in  each  cell.  Deuterium-trit ium  react¬ 
ions  are  counted  during  one  cycle  and,  by  assuming  isotropic 
production  of  one  neutron  per  reaction,  the  new  neutrons  are 
placed  at  the  center  of  the  microsphere  and  attenuated  during 
their  pathlength  through  the  plasma.  Using  an  energy  of 
14.1  Mev  for  the  neutrons  and  a  cross  section  of  0.8  barns 
for  collisions  with  deuterium,  tritium  and  helium  nuclei, 
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the  number  of  interactions  per  cell  is  ascertained.  An 
average  energy  transfer  per  collision  is  used  to  deposit 
energy  into  a  cell  based  on  the  specific  ion  numbers  in 
the  cell. 

5.  Subroutine  ALPHA1 .  This  subroutine  is  designed  to 
compute  the  energy  deposited  in  each  cell  by  monoenergetic 
3.5  Mev  alpha  particles  born  in  the  deuterium-tritium 
reactions.  Additionally,  it  computes  the  number  of  alphas 
produced  in  each  zone,  adjusts  the  helium  particle  mass 
number  in  each  zone,  and  adjusts  ion  and  electron  temperature 
in  each  zone.  Collisions  with  electrons  are  assumed  to 
dominate  the  alpha  particle  deceleration  and  are  also 
assumed  to  provide  no  significant  scatter  to  the  alpha  parti¬ 
cles.  Upon  reaching  an  area  where  ionic  collisions  dominate 
energy  loss,  the  motion  of  the  alpha  particle  is  stopped 
and  all  residual  energy  is  dumped  in  the  current  cell. 

Subprogram  ALPHA1  calls  14  other  subroutines  to 
complete  its  tasks.  General  program  organization  is  seen 
in  Figure  3. 

Subroutine  ALPHA!  updates  the  total  heat  added,  the 
number  of  alpha  particles,  total  mass,  electron  temperature, 
and  ion  temperature  in  each  cell.  Geometry  subroutines 
compute  alpha  particle  position,  direction  and  cell  number 
as  the  particle  slows.  An  integration  increment  is  computed 
to  integrate  dE/ds  along  the  path  length  in  each  cell  using 
a  4th  order  Runge-Kutta  method.  The  energy  deposited  in 
each  cell  is  partitioned  between  ions  and  electrons. 
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Figure  3.  General  Program  Organization  of  Subroutine  ALPHA1 

6.  Subroutine  HTFLX.  This  subroutine  computes  the 
energy  transfer  between  the  radiation,  electron  and  ion 
components.  Drawing  heavily  on  Appendix  B  of  Reference  9, 
HTFLX  first  models  electron-ion  energy  exchange  and  electron 
heat  conduction.  Ionic  heat  conduction  is  considered  negligible. 
Energy  exchange  between  the  photons  and  electrons  is  then 
calculated.  The  energy  exchange  between  photons  and  ions 

is  considered  negligible. 

7.  Subroutine  OUTPUT.  This  subroutine  prints  a  variety 
of  information  on  hard  copy  including  iteration  number,  time. 
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time  increment  used  during  the  iteration,  input  energy,  and 
output  energy.  Additionally,  OUTPUT  provides  current  cell 
radius,  velocity  of  the  cell  wall,  reactions  during  the 
iteration,  relative  density,  electron,  ion,  and  radiation 
temperatures,  energy  output,  heating,  and  neutron  fluence 
for  each  cell.  Subroutine  OUTPUT  can  be  called  at  convenient 
places  in  the  main  program  to  provide  timely  output  of  data. 
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Approach 

This  report  will  briefly  describe  phenomena  important 
to  nuclear  fusion  for  a  deuterium-tritium  fuel.  Energy 
transfer  processes  in  the  plasma  will  then  be  addressed 
including  Coulomb  collisions,  bremsstrahlung  and  Compton 
scatter.  The  equations  used  to  code  MOXNEX  are  then 
discussed  as  each  subroutine  is  examined  in  detail. 

Chapter  IV  discusses  the  validation  of  the  hydro¬ 
dynamics  using  a  quasiequilibrium  study  and  a  point  explo¬ 
sion.  Chapter  V  addresses  recommendations  and  code  improve¬ 
ments  . 

Appendices  include  a  program  listing,  a  glossary  of 
variables,  equation  derivations,  and  discussions  of  pertinent 
items.  Each  is  referred  to  in  the  text  where  appropriate. 
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Introduction 


II. ‘  Physical  Processes 


Before  looking  at  the  MOXNEX  code,  some  basic  phenomena 
concerning  thermonuclear  plasmas  and  fusion  will  be  examined. 
The  question  of  why  fusion  is  attractive  will  be  addressed, 
as  will  the  nuclear  reactions  of  interest.  Knowing  why  the 
process  is  attractive,  attention  will  be  turned  to  the 
question  of  feasibility.  A  simple  three  species  plasma 
model  is  then  presented  and  energy  transfer  processes  between 
these  species  examined.  The  three  species  plasma  model  is  then 
generalized  to  a  thermonuclear  three  species  model  that  describes 
all  MOXNEX  processes. 

Fusion 

As  mentioned  previously,  energy  resources  are  a  growing 
national  concern.  The  following  paragraphs,  then,  will 
reveal  why  energy  can  be  attained  from  a  fusion  technology. 

The  average  mass  per  nucleon  varies  slightly  from  1 
atomic  mass  unit  in  most  nuclei.  Table  1  gives  values  for 
the  average  mass  per  nucleon  of  selected  light  nuclei.  •  If 
nuclear  reactions  can  be  produced  such  that  the  total  mass 
of  the  reactants  is  greater  than  the  total  mass  of  the 
products,  a  net  release  of  energy  of  E  =  (Am)c2  is  attained. 
The  quantity  Am  is  known  as  mass  defect. 
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Table  1 


Average  Mass  per  Nucleon  for  Selected  Light  Nuclei 

(Refs  16  and  5) 


For  instance,  consider  the  reaction 

D+  T-*  H«4  +  r>  cd 

This  will  probably  be  a  central  reaction  in  a  first  generation 
fusion  reactor  (Ref  6:279).  The  total  nuclear  mass  of  deuterium 

t 

and  tritium  is  5.029254604  amu  and  the  total  mass  of  a  helium-4 
nucleus  and  a  neutron  is  5.010171114  amu.  The  difference  is 
0.0188349  amu.  This  equates  to  an  energy  release  of  17.609889 
Mev.  The  process  of  exothermically  combining  nuclei  is  known 
as  fusion. 

Although  many  fusion  reactions  are  .possible,  the  primary 
reactions  of  interest  for  hydrogen  fusion  are  listed  in 
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Table  2.  The  excess  kinetic  energies  of  the  reaction 
products  are  indicated  with  the  reaction  products. 

Table  2 


Reactions  of  Interest  for  Hydrogen  Fusion 

(Ref  19:2) 


In  order  to  bring  about  fusion  reactions,  it  is  neces¬ 
sary  that  the  ionized  nuclei  collide  with  sufficient  energy 
to  overcome  Coulombic  repulsion.  Classically,  the  relative 
energy  required  to  overcome  electrostatic  repulsion  for 
hydrogen  is  0.28  Mev  (Ref  10:7).  The  quantum  mechanical 
effect  of  barrier  penetration,  however,  reduces  this  energy 
threshold. 

The  nuclei  are  distributed  in  kinetic  energy  or  velocity 
in  a  Maxwellian  distribution.  The  average  kinetic  energy 
of  a  particle  in  a  Maxwellian  distribution  is  y  kT^. 

There  are  decreasing,  but  still  finite,  numbers  of  particles 
at  many  kT^,  however,  and  the  high  energy  of  the  particles 
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in  the  "Maxwellian  tail"  is  therefore  a  notable  influence  in 
overcoming  Coulombic  repulsion.  These  two  factors,  barrier 
pentration  and  the  high  energy  tail  on  the  Maxwellian 
distribution  of  the  ions,  make  consideration  of  fusion  reac¬ 
tors  feasible. 

Again,  for  a  feasible  system,  the  fusion  reaction  rate 
must  be  at  least  finite.  Generally,  a  reaction  rate  may 
be  expressed 

RR=n,nxJ-v  (2) 

where 

RR  is  the  reaction  rate  density 
n^  is  the  number  density  of  species  1 
n2  is  the  number  density  of  species  2 
< T  is  the  reaction  cross  section 
v  is  the  relative  velocity  of  the  reactants. 

The  nuclear  cross  sections  in  Eq  (2)  are  strong  functions 
of  velocity  or  kinetic  energy.  Averaging  over  the  entire 
range  of  relative  velocities,  the  product  <TV  can  be 
represented  as  an  expectation  value  .  Values  of  <<r*> 

which  assume  a  Maxwellian  particle  distribution  are  shown 
in  Figure  4.  Notable  from  Figure  4  is  the  deuterium-tritium 
reaction  is  much  more  likely  than  the  deuterium-deuterium 
reaction  at  any  given  kinetic  temperature. 
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The  excess  kinetic  energy  may  be  redeposited  in  the 
plasma  or  it  may  escape  from  the  mixture.  In  order  to 
sustain  the  thermonuclear  burn,  the  redeposition  of  energy 
into  the  plasma  is  preferable.  If  more  energy  can  be  retained 
in  the  plasma,  less  energy  is  required  from  an  outside  input. 
This  energy  may  be  transferred  from  the  reaction  products  by 
Coulomb  collisions,  participation  in  elastic  collisions,  or 
radiation-electron  energy  transfer.  The  process  of  redeposit¬ 
ing  the  particle  energy  released  by  fusion  reactions  back 
into  the  plasma  is  known  as  bootstrap  heating. 

Singling  out  the  deuterium- tritium  reaction, 

D  +  T-*  He*  (3.5  Mav)  +  n(f4 .(  Niev)  C3) 

the  neutrons  can  participate  only  in  elastic  collisions,  but 
the  alpha  particles  can  participate  in  all  three  mechanisms 
named  above. 

The  alpha  particle  path  is  dominated  by  Coulomb 
collisions  with  both  electrons  and  ions.  Thus,  much  or  all 
of  the  3.5  Mev  possessed  by  the  alpha  particle  is  retained 
in  the  plasma  promoting  further  fusion. 

In  order  to  retain  the  14.1  Mev  of  the  neutron  within 
the  plasma,  high  plasma  densities  are  necessary  (Ref  4:330). 

The  mean  free  path  of  a  14.1  Mev  neutron  at  typical  magnetically  con¬ 
fined  thermonuclear  plasma  densities  CIO1 5  ions/cm3)  is  ~109  cm3  so 
neutrons  are  difficult  to  retain.  The  energy  transfer  from 
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radiation  energy  is  also  generally  a  small  factor.  The 
photon  mean  free  path  for  particle  number  densities  10*5  cm^ 
is  ^10"°  cm  (Ref  10:29).  Significant  bootstrap  heating 
comes  primarily  from  alpha  particle  Coulomb  collisions  for 
the  reaction  given  in  Eq  (3).  A  generalization  can  be  made  to 
other  reactions  that  significant  bootstrap  heating  comes 
primarily  from  the  charged  particle  reaction  products. 

A  primary  energy  loss  mechanism  is  bremsstrahlung .  As 
ionized  particles  are  deflected  by  Coulomb  interactions, 
they  continuously  radiate.  The  rate  at  which  energy  is 
radiated  by  an  accelerating  charge  of  z  =  1,  Pr,  is  expressed 
by  the  Larmor  formula,  namely, 

P_  .  2,  e*  <=(■*• 

r  3  c3 

where 

e  is  the  electronic  charge  in  statcoulombs 
a  is  the  acceleration  in  cm/sec2 
c  is  the  speed  of  light  in  cm/sec. 

This  newly  created  photon  energy  can  be  redeposited  through 
inverse  bremsstrahlung  or  Compton  collisions,  but  it  may 
also  escape  the  plasma.  An  optically  thick  plasma  is 
desirable  to  recapture  this  radiation. 

If  the  rate  of  energy  loss  from  the  plasma  due  to 
bremsstrahlung  is  greater  than  the  rate  of  energy  deposition 
by  bootstrap  heating,  energy  must  be  put  into  the  system 
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from  an  outside  source  in  order  to  sustain  the  thermonuclear 
burn.  If  the  energy  deposition  rate  is  greater  than  brems- 
strahlung  loss  rate,  reactions  may  continue  as  long  as  fuel 
is  available.  A  bremsstrahlung  energy  loss  rate  may  be  calcu¬ 
lated  as  a  function  of  number  density  and  temperature.  Postu¬ 
lating  that  100%  of  the  charged  particle  energy  is  retained 
in  the  plasma,  a  bootstrap  heating  rate  may  also  be  calculated 
as  a  function  of  number  density  and  temperature. 

Fortunately,  as  ion  temperature  increases,  the  bootstrap 
heating  rate  increases  faster  than  the  bremsstrahlung  loss 
rate  (Ref  10:33-36).  The  temperature  at  which  the  bootstrap 
heating  rate  equals  the  bremsstrahlung  loss  rate  is  known  as 
minimum  ideal  ignition  temperature.  Below  this  temperature, 
energy  must  be  pumped  into  a  system  to  realize  sustained  thermo 
nuclear  burn.  Above  this  temperature,  the  reaction  may  be 
self-sustaining  and  energy  can  be  extracted  from  the  system. 

Noting  the  charged  particle  products  of  the  reactions 
from  Figure  2  and  the  probability  of  reactions  from  Figure  3, 
this  "ideal  ignition  temperature"  for  deuterium-tritium  is 
about  4  kev  and  for  deuterium-deuterium  is  about  36  kev, 
almost  an  order  of  magnitude  difference  (Ref  10:35)., 

These  temperatures  correspond  to  particles  with  high 
kinetic  energies  which  tend  to  diffuse  the  plasma.  Coulomb 
repulsion  adds  to  dispersion  effects.  Confinement  of  the 
plasma  for  times  long  enough  to  achieve  an  economical  energy 
gain,  therefore,  becomes  difficult.  Since  the  power  density 
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depends  on  the  reactant  number  density,  the  confinement 
time  for  energy  breakeven  is  a  function  of  the  reactant 
number  density.  This  can  be  stated  in  the  Lawson  criteria, 

n'T'  ■S;  lo1*  (5) 

where  T  is  the  confinement  time  (Ref  6:282). 

Magnetic  fusion  relies  on  various  geometries  of  magnetic 
bottles  to  confine  the  ionized  particles.  Typical  number 
densities  for  magnetic  fusion  are  about  10^  particles  per 
cubic  centimeter.  This  requires  a  confinement  time  of  about 
0.1  seconds  to  achieve  energy  breakeven. 

Inertial  confinement,  on  the  other  hand,  realizes 

2  f) 

number  densities  of  about  10  particles  per  cubic  centi¬ 
meter  and  so  requires  a  confinement  time  of  only  about  10 
seconds  (Ref  6:282). 

Additionally,  fusion  reactors  have  some  advantages  over 
fission  reactions.  Deuterium  for  fusion  fuel  is  readily 
available  in  seawater.  Though  tritium  is  not  present  in 
large  enough  quantities  to  be  of  interest  in  seawater  and 
must  be  processed  from  other  sources,  deuterium  is  present 
in  a  ratio  of  about  1  deuterium  atom  to  6500  hydrogen 
atoms  in  ordinary  water  (Ref  10:2). 

Another  advantage  is  that  radioactive  reactants  and 
by-products  are  minimized  using  fusion  when  compared  with 
fission.  Neutrons  do  activate  some  nuclei  in  the  local 
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environment  and  some  materials  desirable  for  fusion  reactions 
are  naturally  radioactive;  for  instance,  tritium.  Fusion 
reactions,  however,  create  no  radioactive  fragments  which 
are  a  necessary  by-product  of  fission. 

In  summary,  these  paragraphs  have  tried  to  briefly 
develop  some  of  the  important  concepts  for  nuclear  fusion  and 
fusion  technology.  Mass  defect,  barrier  penetration,  the 
high  energy  "tail"  of  the  Maxwellian  distribution,  reaction 
rates,  bootstrap  heating,  bremsstr ahlung  losses,  ideal  ignition 
temperature,  and  the  Lawson  criteria  are  all  concepts  central 
to  nuclear  fusion  and  fusion  reactor  technology.  A  quantita¬ 
tive  physical  model  is  now  required  that  can  follow  some  of 
these  concepts  in  a  time  history. 

Three  Temperature  Model 

Knowing  some  of  the  central  concepts  of  fusion,  a  physical 
model  is  now  required  that  can  follow  macroscopic  quantities 
through  a  time  history  of  a  fully  ionized  plasma.  A  simple 
yet  physically  accurate  model  is  desired. 

The  motivating  basis  for  researching  nuclear  fusion  is 
that  it  offers  the  possibility  of  net  energy  gain.  Energy 
may  be  calculated  from  power  over  a  complete  time  history 
and  power  may  be  computed  from  number  densities  and  particle 
velocities  or  particle  energies.  Ionic  hydrogen  will  be  the 
fusion  reactants  in  the  MOXNEX  code  so  ion  number  densities 
and  ion  energies  will  be  required  to  compute  overall  energy 
output.  Ion  number  densities  can  be  characterized  by  an 


energy  distribution.  An  ion  species  will  be  characterized 
by  a  Maxwellian  distribution  which  can  be  denoted  by  a 
number  density  per  unit  energy  at  a  given  kinetic  temperature. 
Ion-ion  collisions  partition  energy  throughout  the  species  to 
attain  a  Maxwellian  distribution  (Ref  26:136). 

But  ion-ion  collisions  are  not  the  only  collisions 
occurring.  Fuel  pellet  atoms  all  bind  one  electron  prior  to 
the  driver  pulse.  Complete  ionization,  then,  means  as  many 
electrons  are  present  as  are  ions.  Electron-electron  col¬ 
lisions  partition  energy  in  the  electron  species  also  in  a 
Maxwellian  distribution  (Ref  26:136).  Electron  number  density 
can  also  be  identified  as  a  function  of  kinetic  temperature. 
This  temperature  may  differ,  however,  from  the  ion  temperature 
(Ref  26:136). 

If  these  temperatures  do  differ,  equilibration  between 
species  will  take  place  through  electron-ion  collisions. 

Energy  exchange  between  these  species  must  also  be  tracked. 

As  mentioned  in  the  previous  section,  Coulombic  accel¬ 
eration  of  the  electrons  by  the  ions  will  result  in  brems- 
strahlung  so  a  photon  species  is  present  also.  A  third 
species  energy  distribution  is  required.  These  photons  may 
be  characterized  by  a  Planckian  distribution  which  requires 
a  third  species  temperature.  Additionally,  the  radiation 
may  couple  energy  back  into  the  electron  species  by  either 
inverse  bremsstrahlung  or  Compton  scatter. 


The  minimum  requirements,  then,  are  two  Maxwellian 
distributions  and  a  Planckian  distribution  characterized 
by  three  species  temperatures.  All  species  can  interact 
with  each  other,  although  the  photons  will  interact  preferentially 
with  the  electrons  because  of  the  mass  difference  between  the 
ions  and  electrons.  In  addition  to  energy  flow  between  species, 
energy  flow  in  space,  diffusion,  will  affect  the  energy  distri¬ 
butions.  The  more  massive  ions  will  move  relatively  slowly 
compared  to  the  lighter  electrons,  but  electron  diffusion  must 
be  considered.  Also,  radiation  diffusion  will  carry  energy 
through  space. 

A  simple  model  is  now  apparent.  A  three  temperature 
model  using  separate  ion,  electron  and  radiation  temperatures 
accounting  for  electron  and  radiation  diffusion  may  be  used 
to  model  a  fully  ionized  plasma.  Such  a  model  is  illustrated 
in  Figure  5, 

The  task  at  hand  now  is  to  develop  understanding  of  the 
pertinent  energy  transfer  processes.  This  will  be  the  topic 
of  the  next  section. 

Coulombic  Phenomena 

Electrostatic  encounters  between  particles  of  charge 
and  Z ^  which  are  separated  by  a  distance  r  are  described 
by  the  Coulomb  potential,  namely 


+  CLC.CTTROUS 


f 


BKCMSSTRAUIWCi  ' 

I  1NVER5C.  i 

©RjCrlSSTl^LUMG 
\  A  £jOM?>TOi4 

\  T  C0Ua£>10*4f> 


COUUO^fb 

cov-i-\sio»ao 


f  I  motion 


IONS 


FlCdJKE  5  :  CM&ZO+  “TRANSFER.  BETTWEEN  SPCCIES 
IN  A  FULIY  ION1ZCP  PLA3MA. 


I 


In  a  plasma,  these  electrostatic  forces  contribute  signi¬ 
ficantly  to  energy  transfer  between  species  and  to  species 
diffusion. 


r 


The  distance  r  is  a  determining  factor  in  the  magnitude 
of  the  resulting  force.  In  the  case  of  charged  particles 
possessing  kinetic  energy,  an  impact  parameter,  b,  is  often 


used  to  denote  the  distance  of  closest  approach.  Figure  6 
illustrates  the  impact  parameter  of  approaching  particles 
in  the  rest  frame  of  particle  2. 


Figure  6.  Illustration  of  Impact  Parameter 
Between  Two  Charged  Particles 
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Figure  5  assumes  a  small  momentum  transfer  so  that  particle 
1  is  basically  undeflected  from  straight  line  motion.  In  a 
center  of  mass  frame,  the  paths  of  both  particles  is  actually 
hyperbolic. 

In  space,  the  range  of  the  Coulomb  force  is  infinite  and 
there  is  no  maximum  interaction  distance.  In  a  plasma, 
however,  even  though  the  entire  system  is  electrically  neutral, 
random  thermal  motion  creates  small  deviations  from  strict 
electrical  neutrality.  Electrostatic  forces  accentuate  this 
and  a  time  average  will  show  a  net  negative  charge  density 
around  an  ion  and  a  net  positive  charge  density  around  an 
electron.  This  net  opposite  charge  may  be  regarded  as  an 
"atmosphere”  around  each  charged  particle  (Ref  10:84).  The 
radius  of  this  atmosphere  based  on  a  Maxwellian  distribution 
and  a  Coulomb  potential  expressed  in  Poisson's  equation  is 


(Ref  6:10).  The  quantity  ^  is  known  as  the  Debye  length 
and  determines  the  maximum  shielding  or  screening  distance 
at  which  a  Coulomb  force  can  be  felt  in  a  plasma.  Notice 
that  as  electron  temperature  increases,  the  Debye  length 
increases.  It  is  also  a  function  of  electron  number  density. 
The  Debye  length  is  a  dynamic  parameter  and  the  pseudo 
"atmosphere"  is  not  static. 
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Coulomb  collisions  can  be  divided  into  short  range 
encounters  and  long  range  encounters.  The  Debye  length  is 
the  limiting  distance  for  long  range  encounters.  The 
charged  particle  is  "screened"  from  charged  particles  farther 
away  because  of  the  opposite  charge  "atmosphere."  Because 
of  the  large  number  of  particles  within  a  Debye  radius,  the 
number  of  long  range  encounters  is  large.  The  number  of 
long  range  encounters  is  so  great  that  their  effect  greatly 
outweighs  the  effect  of  short  range  collisions  (Ref  26:123). 

The  cross  section  for  a  short  range  encounter  resulting 
in  an  angle  of  deflection  from  its  original  path  of  90° 
is 
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where  E-^  is  the  energy  of  the  less  massive  particle  in  the 
rest  frame  of  the  more  massive  particle.  In  a  hydrogen 
plasma,  this  can  be  an  electron  passing  any  positive  particle. 

A  long  range  encounter  may  be  described  by  a  change  in 
momentum  such  that 
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This  is  equivalent  to  multiple  scatters  summing  to  a  total  90* 
deflection.  The  cross  section  for  this  type  of  interaction  is 
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where  is  again  the  energy  of  the  less  massive  particle 

in  the  rest  frame  of  the  more  massive  particle,  and 
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or 
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Both  short  range  and  long  range  Coulomb  cross  sections  are 
discussed  in  Appendix  C. 

Comparing  Eqs  (8)  and  (10),  note  that  the  long  range 
cross  section  is  =  8£nA  greater  than  the  short  range  cross 
section.  Also  noting  that  typical  values  for  tnA  are  between 
10  and  25  at  thermonuclear  temperatures  and  densities 
(Ref  10:94),  the  long  range  interaction  is  much  more  dominant. 

Relaxation  times  are  often  used  to  describe  long  range 
encounters  and  may  be  defined  as  the  time  required  for  the 
momentum  changes  to  sum  such  that 
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The  long  range  cross  section  given  in  Eq  (10)  yields  for 


electron-ion  relaxation  time 


4  TT  eA  r\i  /sr\ 


(14) 


Coulomb  interactions  between  ions  and  electrons  at 
different  temperatures  will  eventually  lead  to  equilibrium. 
Based  on  Maxwellian  velocity  distributions  for  both  ions 
and  electrons  and  defining  t  from  the  relation  given  by 
Spitzer  (Ref  26:135), 
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where  T-  is  species  kinetic  temperature,  t  can  be 
J 

written 

t  =  3m.  j  T]_ 
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Relaxation  times  and  equilibration  times  are  further 
discussed  in  Appendix  N. 

In  addition  to  electron- ion  energy  exchange,  Coulomb 
collisions  add  to  energy  transfer  through  heat  flow.  Heat 
conduction  may  be  approached  in  the  same  manner  as  other 
transport  phenomena  using  a  conductivity  coefficient 
(Ref  26:143). 

In  the  presence  of  a  temperature  gradient  (assuming 
no  external  electric  field),  the  flow  of  heat,  Q,  can  be 
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described  by 


Q  -  KVT 


cm 


where  T  is  temperature  and  K  is  the  conductivity  coefficient.  For 
a  Lorentz  gas,  a  fully  ionized  gas  in  which  the  electrons  are  assumed  to 
interact  only  with  ions  which  are  all  at  rest  (Ref  26:138),  this  coeffi¬ 
cient  is  given  by  Spitzer  (Ref  26:144), 
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Heat  flow  caused  by  a  temperature  gradient  results  in  a 
current,  however,  and  this  current  produces  a  secondary 
electric  field  which  reduces  the  flow  of  heat  by  a  factor  of 
e  which  is  dependent  on  the  net  Z  of  the  material. 
Additionally,  for  an  actual  gas  K  is  further  reduced  by  a 
factor  of  6t  which  is  also  dependent  on  the  net  Z  . 
Employing  these  terms,  the  conductivity  coefficient  becomes 
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A  conductivity  coefficient  may  be  derived  for  any  charged 
species.  Equation  (19)  gives  this  coefficient  for  the 
electron  species.  As  the  least  massive  charged  particles  in 
a  thermonuclear  plasma,  the  electrons  should  dominate  heat 
transport  if  the  electron  temperature  -is  not  very  much  smaller 
than  the  ion  temperature. 
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Bremsstrahlung  Processes 

Bremsstrahlung  occurs  when  a  particle  having  charge  and 
finite  kinetic  energy  is  accelerated,  resulting  in  a  photon 
being  radiated  and  decreased  kinetic  energy  of  the  original 
moving  particle.  This  occurs  typically  in  a  thermonuclear 
plasma  as  an  electron  is  accelerated  as  it  passes  a  positive 
ion.  The  instantaneous  power  radiated  in  a  non-relativistic 
acceleration  is  given  by  the  Larmor  formula  (Ref  14:469), 
here  for  Z  =  1, 
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A  bremsstrahlung  emission  is  shown  symbolically  in  Figure  7. 
for  the  case  of  an  electron  passing  a  positive  ion. 


Figure  7.  Symbolic  Illustration  of  Bremsstrahlung 


Kith  these  pictures  in  mind,  some  of  the  theory  pertinent 
to  bremsstrahlung  interactions  will  now  be  developed.  The 
equations  presented  in  the  following  discussion  and  Appendices 
E  through  G  are  prompted  by  the  unpublished  notes  of  Dr. 

George  H.  Nickel.  Significant  portions  of  the  development 
are  contained  in  Reference  22. 

The  bremsstrahlung  cross  section  for  a  non-relativistic 
Coulomb  collision  is  given  by  Jackson  (Ref  14:513)  as 


where 

is  the  charge  of  the  radiating  particle 
Z2  is  the  charge  of  the  particle  providing  the 
accelerating  force 

M  is  the  mass  of  the  radiating  particle 
v  is  the  velocity  of  the  radiating  particle 
or  E  in  the  rest  frame  of  the  ion. 

This  cross  section  is  also  supported  by  Bethe  and  Heitler 
(Ref  14:512).  Its  dimensions  are  area  per  energy. 

For  the  case  of  a  thermonuclear  plasma,  the  radiating 
particle  is  usually  an  electron  so  that  Z  =  1  and  M  =  me 
and  the  particle  providing  the  accelerating  for  ionic  hydrogen 
or  ionic  helium.  Referring  to  Appendix  E  and  denoting  the 
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ionic  charge  as  simple 


Z  ,  Eq'(20)  becomes 
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where 


ok  =  e'/'fic  is  the  fine  structure  constant 
Eg  is  the  electron  energy 
E  is  the  photon  energy 


The  bremsstrahlung  reaction  rate  density  is  then 


RR*  -  ntr\t(E*)'ir rt(Ee,E i) 


(22) 


which  has  dimensions  of  reactions  per  volume  per  time  per 
square  energy.  The  inverse  bremsstrahlung  cross  section  may 
be  derived  using  this  reaction  rate  density  and  assuming 
local  thermodynamic  equilibrium  exists.  The  inverse  brems¬ 
strahlung  cross  section,  also  discussed  in  Appendix  E,  is 
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with  dimensions  area  per  energy. 

The  reaction  rate  densities  for  the  respective  processes 
can  then  be  expressed  using  these  cross  sections.  The  brems¬ 
strahlung  reaction  rate  density  can  be  written  explicitly  as 
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which  has  dimensions  of  reactions  per  volume  per  time.  From 
this,  the  power  density  can  be  expressed  as 
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which  has  dimensions  of  energy  per  unit  volume  per  unit  time. 
The  enhancement  factor 
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used  for  the  bremsstrahlung  equations  is  derived  in  Appendix  D. 

For  inverse  bremsstrahlung,  the  reaction  rate  density 
and  the  power  density  can  be  expressed  as 

'Op  r  co 
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Substituting  explicit  expressions  into  ,  changing 

variables  to 
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(29) 


(30) 


and  noting 


(31) 


results  in 


(32) 


Analogous  manipulations  can  be  performed  to  refine 
This  is 

e*>  I  \  r0D 
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These  two  power  densities  can  be  expressed  using  a 
Gaunt  factor  G C JO ^  »  where  j  denotes  either  pure  brems- 
strahlung  or  pure  inverse  bremsstrahlung .  Using  this  Gaunt 
factor,  Eqs  (32)  and  (33)  may  be  written  as 
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(34) 


and 


Recognize  that  these  power  densities  are  limiting  cases. 
Equation  (34)  models  a  case  of  pure  bremsstrahlung  and 
Eq  (36)  models  the  case  of  pure  inverse  bremsstrahlung.  For 
practical  applications,  the  power  density  equations  must  be 
extended  from  the  limiting  cases  of  either  pure  bremsstrahlung 
or  pure  inverse  bremsstrahlung  to  a  situation  of  competing 
processes.  This  is  done  by  simple  subtraction. 

Note  that  the  power  density  equations  have  a  common 
coefficient,  which  may  be  denoted  by  ,  and  is 
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Substituting  explicit  expressions  for  e<  ,  rQ  ,  ng  ,  and 


n.  ,  this  becomes 
1  ’ 


A*  =  32/2  \'/a  e4N,  ,  (?'\  6*  £J_  r(v) 
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Thus  the  net  rate  of  energy  flow  into  the  radiation  field 


from  bremsstrahlung  processes  is 
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Denoting 
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this  can  be  written 
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where 


f(0=  C^(v/7«^)eWxcJ* 

'-'i 

Note  that  G()f)  is  dimensionless. 

The  values  of  G()f)  as  ^  0  and  as  oc>  may  be 

found  analytically.  These  are 

G(o)  =  1 


G(eo)  = 


Appendix  F  discusses  the  reaction  rate  densities,  power 
densities.  Gaunt  factors,  and  the  bremsstrahlung  coupling 
coefficient  in  more  depth.  Appendix  G  discusses  the  Gaunt 
factor  and  its  limiting  values. 
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The  quantum  theory  of  light  postulates  that  photons 
behave  like  particles  except  for  the  absence  of  rest  mass. 
This  concept  provides  a  foundation  for  radiation-material 
energy  interaction  to  be  treated  as  an  elementary  mechanical 
collision. .  Figure  9  symbolically  illustrates  such  a 
collision  known  as  a  Compton  collision. 


Figure  9  .  Symbolic  Illustration  of  a  Compton  Collision 


The  energy  transfer  can  be  described  by  the  wavelength  change 
of  the  photon  which  is  (Ref  2:71) 


Since  the  photon  is  not  absorbed  and  a  photon  exits  from 
the  interaction,  Compton  collisions  are  a  scattering  process 
as  opposed  to  absorption  processes  such  as  pair  production 
or  photoelectric  effect. 

The  classical  cross  section  for  Compton  scatter  is 

crT-  — 


where 


r.  .-Sl 


is  the  classical  electron  radius.  This  cross  section,  known 
as  the  Thomson  cross  section,  becomes  suspect  at  higher 
photon  energies. 

At  higher  photon  energies,  Compton  scatter  can  be  des¬ 
cribed  by 

^ =  4'  ■  *6fo) +  t  Ofe)*4  •  •  •] 


(Ref  3:124),  where  Ev  is  the  energy  of  the  incident  photon. 
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This  is  known  as  the  Klein-Nishina  cross  section  and  applies 
only  for  photon  energies  such  that  mg  c^ 

Inverse  Compton  collisions  are  also  possible  and  such  an 
interaction  is  symbolically  represented  in  Figure  10. 


Figure  10.  Symbolic  Illustration  of  an  Inverse  Compton  Collision 


Note  that  in  both  processes  a  photon  is  emitted.  Again, 
since  photons  are  Bose  particles,  the  photon  population  will 
be  enhanced  by  photons  already  existing  in  the  plasma  and 
this  enhancement  factor  can  be  expressed  as 
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Again,  much  of  the  following  development  is  due  to 
George  H.  Nickel  (Ref  22). 

The  reaction  rate  density  for  Compton  scatter  now  can 
be  written  as 


r>y  (Ev)  a; 


(53) 


which  has  dimensions  of  reactions  per  volume  per  time  per 
energy.  This  formula  assumes  the  electrons  are  at  rest.  The 
cross  section,  however,  is  dependent  on  electron  energy.  Note 
from  Eq  (52)  that  if  Ev«<  mec2  that  the  Klein-Nishina 
cross  section  reduces  to  the  Thomson  cross  section.  This  is 
well  approximated  in  a  thermonuclear  plasma  as  typical  photon 
temperatures  can  be  expected  to  be  on  the  order  of  1  kev 
(Ref  9:475).  As  the  electron  rest  mass  is  mgc2  =  511  kev  , 
an  approximation  might  be  made  assuming  for  a 

thermonuclear  plasma  and  using  for  the  reaction  rate 


A  power  density  may  be  written  from  the  reaction  rate 
by  including  the  energy  exchange  per  collision  and  integrating 
over  the  entire  range  of  photon  energies.  This  average 
energy  exchange  per  collision  can  be  approximated  by 

<^£>  =  I 
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This  quantity  is  discussed  in  Appendix  H. 

The  power  going  from  radiation  to  electrons  per  unit 
volume  is 
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No  integration  is  done  over  electron  energy  as  the  electrons 
are  again  assumed  to  be  at  rest.  Using  explicit  expressions, 


this  becomes 
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By  multiplying  and  dividing  by  (kTr)5  ,  and  substituting 
y  =  E  /k  Tr  , 
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The  value  of  the  integral  is  4^/15  ,  therefore 
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This  is  the  rate  of  energy  flow  from  radiation  to  electrons 
for  pure  Compton  scatter. 

The  calculation  for  inverse  Compton  scatter  is  fairly 
complicated  and  will  not  be  included  in  this  report.  Inverse 
Compton  requires  that  the  electron  have  non-zero  kinetic 
energy  prior  to  incidence  with  the  photon.  The  calculation 
first  requires  a  transformation  to  the  rest  frame  of  the 
electron.  The  energies  of  the  incident  and  emerging  photons 
are  also  transformed  in  the  electron  rest  frame.  Rates  are 
calculated  in  this  reference  frame  and  transformation  back 
to  the  original  frame  is  required.  The  result  is  (Ref  20) 


This  equation  represents  pure  inverse  Compton  scatter. 

The  net  flow  of  energy  from  electrons  to  radiation  can 
be  found  by  subtracting  the  power  from  radiation  to  electrons 
due  to  Compton  scatter  from  the  power  from  electrons  to  radia 
tion  due  to  inverse  Compton  scatter,  or 
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Inserting  Eqs  (58)  and  (59)  for  P^  and  P^ 

?r-  ^  .... 

The  coupling  coefficient  for  Compton  scatter  may  then  be 
written 


This  can  be  expressed  in  terms  of  the  Stephan-Boltzmann 
constant  and  explicit  expressions  as 
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Pertinent  information  on  the  Compton  power  density  and 
Compton  coupling  coefficient  may  be  found  in  Appendix  I. 

Thermonuclear  Three  Temperature  Model 

Complete  ionization  does  occur  in  a  thermonuclear 
plasma.  Nuclear  interactions,  however,  are  not  modelled  in 
Figure  4.  Neutral  particles  in  the  form  of  neutrons  exist 
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in  a  thermonuclear  state  and  are  also  not  modelled  in 
Figure  5.  Figure  11  includes  these  items. 

The  neutrals  are  a  result  of  the  nuclear  reaction 

D  +  T  -+•  /-/c4(3.saw)  +  r,(t+.t  M*v)  (J) 

Note  that  one  particle  is  returned  to  the  ion  distribution  with 

i 

an  additional  3.5  Mev  of  energy.  One  neutron  at  14.1  Mev  is 
added  to  the  neutral  distribution  for  each  nuclear  reaction, 
also.  The  nuclear  reactions,  then,  affect  both  the  ion  and 
neutral  particle  distributions,  and  are  shown  affecting  both 
in  Figure  11. 

Once  born,  however,  the  14.1  Mev  neutron  may  undergo 
elastic  collisions  with  the  ions,  adding  energy  to  the  ion 
distribution.  This  interaction  is  illustrated  in  Figure  11 
also.  The  neutrons  have  no  mechanism  with  which  to  interact 
with  the  photons  and,  though  the  neutrons  can  collide  elast¬ 
ically  with  the  electrons,  the  collision  cross  section  is 
negligible  when  compared  to  the  ion  collision  cross  section. 

The  energy  distributions  of  the  ions,  electrons  and  photons 
remain  the  same  as  do  the  interactions  between  these  species. 

The  energy  distribution  of  the  neutrons  is  represented  as  a 
Maxwellian  distribution  in  Figure  11.  At  extremely  high 
densities,  this  should  be  the  case  as  neutrons,  being  particles, 
will  tend  toward  a  Maxwellian  distribution.  But  at  lesser 
densities,  many  14.1  Mev  neutrons  will  escape  the  plasma  after 
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only  few,  if  any,  collisions.  Few  will  thermalize.  The 
escaping  particles  will  not  be  in  the  distribution  long, 
however,  so  the  actual  distribution  may  be  a  Maxwellian 
distribution  with  a  peak  at  high  temperatures. 

The  nuclear  interaction  will  dramatically  add  to  the 
energy  of  the  total  system  and  is  indeed  the  motivating 
process  for  this  study.  Chapter  II,  then,  addresses  the 
system  illustrated  in  Figure  11  and  how  these  processes  are 
modelled  in  the  MOXNEX  code. 
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III.  Theory  and  Coding 


Introduction 

This  chapter  discusses  specific  theory  and  how  it  is 
coded  into  program  MOXNEX.  Remembering  energy  transfer 
theory  developed  in  Chapter  II,  the  specific  equations  of 
MOXNEX  will  first  be  presented.  A  comparison  of  this  set 
of  equations  with  another  formalism  will  then  be  given. 

Finally,  the  coding  of  the  formalism  into  the  specific 
equations  of  MOXNEX  will  be  discussed. 

MOXNEX  Formalism 

The  MOXNEX  one-dimensional,  three  temperature,  spherical, 
Lagrangian  code  uses  the  momentum  equation 

_ _ L 

H  f  *r  ^  T  '64> 

where  P^,  is  the  sum  of  the  three  species  pressures  or 

PT  "  pe  +  Pi  +  pr  (65) 


and  q  is  the  artificial  viscosity.  The  coupled  energy 
equations  are 
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The  coupling  coefficients  in  the  equations  are 


a- r  ,  A.k 
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where  j  and  k  denote  specific  species.  Each  energy 
equation  is  derived  using  the  chain  rule 


cK  .  Cv  ii.  + 
6*  ~  VJ  d* 


j  denoting  species,  where 
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and  for  an  ideal  gas 
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The  energy  densities  and  specific  heats  for  the  three  species 
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for  electrons. 
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for  ions, 
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Ion  and  electron  species  pressures  are  computed  using  an 
ideal  gas  equation  of  state.  Radiation  pressure  is 


(80) 


The  ion  electron  fields  are  coupled  by 


where 

V  8  (zrr),/z 

Vc\  3  Az(kTt)** 


(81) 


(82) 


The  electron  diffusion  coefficient  is 
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K  -  20(*  £Jl 

Ke  W)  n?e*ZJr,Mi 


where 


€  Sj  s 


_ 0.432 _ 

(3.44  *  2  +  0.26^121) 


3  \Vl 


—  ■■  v  y, 

(ze>\l  ir N*r J  rz^^c)(m) 


and  is  defined  to  have  a  minimum  value  of  1.  The  electron 
radiation  coupling  coefficient,  Agr  ,  is 

A€r  =  Acr  +  Aer 


where 


a  i9  _  32.  /  z  \/le*Nl_ 

"  3  (ir«J  ^  C 


At t  =  ir 


with 


fiM*I 


«  >r  -K^-dn 

"*{  Until  1-e  J 
L  (!-/)(!-  el/r) 


where 


The  Rosseland  mean  free  path  is 


Ar  = 

f 


from  which  the  radiation  diffusion  coefficient  is  attained, 


namely 


Kr  *  *rTr3 


Only  the  deuterium-tritium  nuclear  reaction  is  modelled 
assuming  a  fuel  that  is  50%  deuterium  and  50%  tritium  so  the 
nuclear  reaction  rate  density  is  modelled  as 


RR, 


=  Ils.  <<rir>B 
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Two  hydrogens  are  subtracted  from  the  hydrogen  number  array 
for  each  reaction. 

Radiation-Electron  Energy  Transfer  Formalism  Comparison 


Reference  9  presents  the  bremsstrahlung  Gaunt  factor 


as 


wy\  f  tft({)Ll-e  J 


(93) 


where 


(94) 


Appendix  F  derives  these  quantities  as 

mi  _  rw^r.-^j 
(<-«( .  -  »•») 


(88) 


where 


f  (0  =  j  A  (VT  4  yjrr)e 


(89) 


All  quantities  are  dimensionless.  The  bremsstrahlung  Gaunt 
factor  is  used  in  the  project  as  derived. 


Reference  9  gives  the  bremsstrahlung  coupling 


coefficient  as 


A 


er 


_  32 /z  Y*e4N*  | 

3  VTt'mJ  fie 


(95) 


and  the  Compton  coupling  coefficient  as 


A 


C 

er 


iii  TreV* 

3  (l*i«  c*)1 


These  are  used  according  to  the  relation 


4T, 

4* 


+ 


where 

A«r  =  4  Act) 


(96) 


(97) 


(98) 


and  j  denotes  either  electron  or  radiation.  The  quantity 

A  „  ,  then,  should  have  dimensions  of  1/mt  .  The  dimensions 
er  ’  ’ 

of  A  in  Eq  (98),  however,  are  l/(t-m-H). 

The  coupling  coefficients  used  in  this  study >may  be 


expressed 


’  ^■v,rS)-N-(ii)r-f‘(T'T 

where 

Aer  *  Aer  +  Aer 


(87) 


(85) 


For  computations  in  MOXNEX 

£H  =:  ±  (Te-To) 

dt  £yj 

where 


(99) 


(100) 


Comparing  Eqs  (86)  and  (95),  it  is  seen  that  they 
differ  by  a  factor  of  prQ  .  Equations  (87)  and  (96) 
differ  by  the  same  factor.  The  dimensions  of  A^r  from 
Eq  (100)  should  be  1/mt  also.  The  dimensions  of  Aer 
in  Eq  (85)  are  l/(t-£s)  .  The  quantity  A^r  ,  then,  is 
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1/mt  ,  consistent  with  Eq  (99)  .  Derivations  of  these 
coupling  coefficients  are  included  in  Appendices  F  and  I. 

Subroutine  GDATA 

This  subroutine  sets  initial  conditions  by  initializing 
variables  and  is  called  one  time  prior  to  the  execution  of 
physical  processes.  Constants  are  gathered  at  the  front  of 
the  subroutine,  then  specific  zone  values  are  computed.  The 
radii  of  the  zones  are  the  inside  radius,  that  is  r^  =  0  . 
This  requires  one  more  radius  element  than  the  number  of 
zones.  Thus,  two  loops  are  necessary;  one  having  an  upper 
limit  of  the  number  of  zones  and  another  to  the  number  of 
zones  +  1. 

Hydrogen  density  and  pusher  mass  are  constants  of  inter¬ 
est  for  parameter  study.  Additionally,  zone  thickness  and 
number  of  zones  may  be  changed. 

Energy  is  introduced  in  the  microsphere  by  initializing 
the  ion,  electron,  and  radiation  temperatures  as  desired. 

If  a  cell  has  no  other  initial  energy  specified,  it  is  set  at 
a  value  of  1.16  °K  or  1  x  10  kev.  Typically,  the  radiation 
temperature  is  set  at  1  kev  and  the  electron  and  ion  temper¬ 
atures  are  set  between  1.8  and  100  kev  following  procedures 
outlined  in  Reference  9. 
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Subroutine  HYDRO 


Subroutine  HYDRO  is  a  one  dimensional  spherical 
Lagrangian  hydrodynamics  code.  It  is  based  on  the  Lagrangian 
energy  equation 


(101) 


where 


T  is  species  temperature 
is  the  specific  heat 
P  is  species  pressure 

q  is  artificial  viscosity 

I  is  species  internal  energy 

V  is  volume 

t  is  time. 


Appendix  J  discusses  this  equation.  Viscosity  is  assumed  to  be  due  only 
to  the  ions,  and  is  modelled  by  artificial  viscosity.  Thus,  q  is  seen 
only  in  the  ion  species  equation. 

The  first  lines  of  the  subroutine  initializes  constants 
and  steps  the  cycle  counter.  Initial  densities  are 

then  calculated  and  viscosities  are  set  to  zero  for  the 
first  cycle.  Ion  and  electron  pressures  are  computed  using 
a  perfect  gas  law  and  species  populations  computed  in  other 
subroutines.  Radiation  pressure  is  calculated  from 
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(80) 


?r-(±£) 


<r  \  T* 


An  effective  temperature  is  calculated  for  the  case 
of  electron  degeneracy  and  is  used  as  a  minimum  electron 
temperature  according  to  the  relation 


tiffed, =  O.COSLS^ 


(102) 


This  is  discussed  in  Appendix  K.  At  low  values  of  Tg  , 
this  effective  minimum  provides  electron  pressures  which 
create  a  disassembly  rate  higher  than  classically  antici¬ 
pated  (Ref  9:477).  Total  pressure  is  the  sum  of  the  species 
pressures,  or 


Pr  =  +  Pe  +  Pr 


The  time  step  is  computed  using  a  fractional  value  of 
the  time  required  for  a  shock  wave  to  travel  across  the 
thinnest  cell  or 


JLt  ..  sr  F 

ttZCAttCfl  * 


(103) 


where 


F  is  the  fractional  value  specified 
c  is  the  speed  of  sound  in  the  medium. 
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The  speed  of  sound  is  computed  using 


This  assumes  a  perfect  gas  so  that  ¥  =  5/3  and  also  that 
viscous  pressure  is  due  only  to  the  ions.  A  minimum  time 
step  is  specified  in  subroutine  GDATA  and  a  maximum  time 
step  is  specified  early  in  the  HYDRO  subprogram.  The 
fractional  value  used  may  be  varied  between  1/10  and  1/2. 

Cell  wall  accelerations  and  velocities  are  computed 
using  total  pressure  and  artificial  viscosities.  New  cell 
wall  positions  are  then  determined. 

With  these  new  positions,  densities  are  updated. 
Artificial  viscosities  are  now  updated  according  to  the 
relation 


<2  -  f>(A  rf(v-v)x 


(105) 


where 

q  is  the  artificial  viscosity 

V«V  is  the  divergence  of  velocity  (Ref  27:136). 

In  the  instance  of  cell  expansion,  artificial  viscosity  is 
set  to  zero. 

Species  temperatures  are  then  updated  accounting  for 
compression  or  expansion.  Computations  assume  the  total 
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Pdv  work  can  be  divided  between  the  ion,  electron 
and  radiation  energies  based  on  the  species  pressures. 

Ion  and  electron  calculations  are  done  using  the  ideal  gas 
law.  Viscosity  is  not  used  in  the  electron  temperature 
update  equation,  again  reflecting  the  assumption  that 
viscous  pressure  is  due  only  to  ions.  Radiation  tempera¬ 
ture  is  updated  using 


-p.  (n+i) 

•  r 


(v^-^sY 


n  'A 


V 


lh+l) 


(106) 


where 

V (n+l)  ^5  tjie  updated  volume 

£V  is  the  change  in  volume 
n  is  the  iteration. 

This  assumes  adiabatic  compression  and  is  discussed  in 
Appendix  M.  Minimum  species  energies  are  specified. 
Finally,  the  total  time  is  re-calculated. 

Subroutine  TBURN 

Subroutine  TBURN  provides  the  reaction  rates  of  the 
deuterium-tritium  fusion  reaction.  Deuterium-deuterium 
reactions  are  not  computed.  Tritium  created  by  deuterium 
deuterium  reactions  is  not  considered.  Constant  tempera¬ 
tures  and  number  densities  of  deuterium  and  tritium  are 
assumed  at  the  time  of  the  subroutine  call. 
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A  minimum  ignition  temperature  flag,  TSTART,  is  set 
at  1  kev,  a  figure  based  on  bremsstrahlung  power  lost.  If 
the  ion  temperature  is  below  the  TSTART  temperature,  no 
reactions  are  computed  in  any  cell.  The  flag  TBURN  is  used 
to  carry  this  information  to  other  subroutines;  TBURN  =  1 
signals  ignition  and  TBURN  =  0  denotes  no  thermonuclear 
burn.  The  number  of  deuterium  and  tritium  ions  burned  are 
subtracted  from  the  totals  in  each  cell. 

The  reaction  of  interest  in  TBURN  is 

P  +  T  ->  He4  +  rv  (i) 


The  reaction  rate  density  for  this  reaction  is 

KRn  '  n>A' do?) 


where  <av>^j  is  the  Maxwell  Boltzmann  velocity  distribution 
weighted  average  of  cross  section  for  the  deuterium-tr itium 
reaction  times  the  velocity  of  approach  of  the  reactants 
(Ref  4:358).  The  weighted  average  of  the  product  <ov>^ 
is  modelled  in  two  ranges,  one  for  ion  temperatures  less 
than  or  equal  to  10  kev,  namely 


£10 


(108) 
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and  for  ion  temperatures  greater  than  10  kev 


<crv) 


/  -i4\  -*h  r(3.638 f/’j-faz"*2'’)] 

"ItohxZ  (3AIX'°  ) ‘  e  (109) 


Values  of  <(<rv>  from  these  functions  are  seen  in 
Figure  11. 

Since  the  original  hydrogen  mix  is  specified  as  50% 
deuterium  and  50%  tritium,  only  the  hydrogen  number  density 
is  used  modelling  the  reaction  rate  density  as 


"RKdt  =  —  <«■*)! 


PT 


(92) 


Subroutine  NUHEAT 

Subroutine  NUHEAT  calculates  total  neutron  heating  in 
each  cell.  The  uncharged  neutrons  do  not  participate  in 
Coulomb  collisions  as  do  alpha  particles  though  neutrons 
can  add  significantly  to  plasma  heating  at  high  fuel 
densities  (Ref  4:330).  Most  neutrons  born  at  14.1  Mev, 
however,  escape  unscattered  for  conditions  of  interest 
(Ref  4:359).  Expense  is  minimized  by  using  a  beam  attenu¬ 
ation  model.  This  simplification  is  an  economical  consider¬ 
ation  balanced  by  the  relative  importance  of  neutron  heating. 
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Figure  12.  Graphical  Values  of  ^crv^vi  Used  in  Subroutine  TRIJRN. 
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Deuterium- tritium  reactions  are  counted  during  one 
cycle  and,  by  assuming  isotropic  production  of  one  neutron  per 


reaction,  the  new  neutrons  are  placed  at  the  center  of  the 
microsphere  and  attenuated  during  their  path  length  through 
the  plasma.  Using  an  energy  of  14.1  Mev  for  the  neutrons 
and  a  cross  section  of  0.8  barns  for  collisions  with  deuter¬ 
ium,  tritium  and  helium  nuclei,  the  number  of  interactions 
per  cell  is  ascertained  according  to 

N HMtfR  or  JwTffcACTlOWS  «■  (TOTAL  NgCnKOH  ^ 

0  “  a  Neutron  "Pehgt Khres.  thi  Ccll) 

The  probability  the  neutron  penetrates  the  cell  is 

p  -(MuHinet.  or  Ions  *  &  Ram/u  "Distance  Throi/sw 

rt‘CN£-rftA1G  *  6 

and  the  ion  number  denisty  is 

The  cross  sections  from  BNL  325  for  deuterium,  tritium 
and  helium  at  14.1  Mev  are  shown  in  Table  3.  Noting  the 
values  in  Table  3,  0.8  barns  seems  low  for  this  treatment. 

The  average  energy  per  collision  is  multiplied  by  the 
number  of  interactions  to  attain  total  heating.  The 
average  energy  per  collision  is  given  by 


(110) 


(111) 


(112) 


Table  3 


Cross  Sections  for  14.1  Mev 
Neutron  Interactions  with  Selected  Nuclei 


nH  +  nHe 


The  factors  of  ^  in  front  of  the  hydrogen  number  densities 
signify  a  50%  mix  of  deuterium  with  a  50%  mix  of  tritium. 

Finally,  subprogram  NUHEAT  calculates  the  number  of 
3.5  Mev  alpha  particles  that  will  give  an  equivalent  amount 
of  heating  if  all  3.5  Mev  is  deposited  in  the  plasma.  This 
number  density  of  alpha  particles  is  passed  in  the  NUHEAT 
argument  list  to  other  subroutines. 

Subroutine  ALPHA1 

The  subroutine  is  designed  to  compute  the  energy 
deposited  in  each  cell  by  monoenergetic  3.5  Mev  alpha 
particles  born  in  deuterium-tritium  reactions.  Also,  it 
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computes  the  number  of  alpha  particles  produced  in  each 
zone,  recomputes  helium  particle  number  density,  recalculates 
total  zone  mass  and  adjusts  ion  and  electron  temperatures. 

Subprogram  ALPHA1  calls  14  other  subroutines  to  accom¬ 
plish  tasks.  Subprogram  organization  is  illustrated  in 
Figure  13. 


Figure  13.  Subprogram  ALPHA1  Substructure 


The  main  subroutine,  ALPHA1,  is  the  location  where 
temperatures,  zone  mass,  alpha  number  densities,  and  heating 
are  computed.  Subroutines  XYT3  through  XYT9  and  GETA  through 
GETB  are  geometry  subroutines.  Subroutine  FINDDS  calculates 
an  integration  increment  and  ELOSS  computes  energy  deposition 
and  energy  partition.  A  logical  flag,  IALPH,  is  set  according 
to  remaining  alpha  particle  energy  and  signals  particle 
behavior.  Table  4  defines  parameter  IALPH.  Before  each  call 
IALPH  is  set  to  0.  This  value  is  changed  by  the  geometry 
subroutines.  The  logical  flag  value  ISTOP  =  0  signals  thermal 
ization. 


Table  4 

Logical  Flag  IALPH  Definitions 


IALPH  Value  Behavioral  Treatment  of  Alpha  Particle 


2  Particle  travels  line  of  sight  until 

thermalization;  electron  collisions 
dominate 

1  Particle  travels  line  of  sight  until 

ion  collisions  dominate;  scatter 
becomes  significant  and  particle  stops 


0 


Particle  stops  in  birth  zone 


Figure  14.  Angular  Directions  for  Alpha  Particle 

1"  Vi 

Leaving  the  j  Zone 

The  geometry  used  by  the  coding  is  summarized  in 
Figure  14.  Alpha  particles  are  assumed  to  be  at  zone 
center  at  birth.  Angles  are  most  finely  divided  at  ir/2  . 
The  paths  are  weighted  by  the  amount  of  angle  considered 
assuming  isotropic  propagation  from  particle  birth  location. 
The  purpose  of  the  geometry  subroutines  is  to  find  the 
distance  traveled  within  a  zone,  thus  enabling  energy 
deposition  modelling.  The  particle  direction  is  assumed  to 
be  unchanged  by  collisions  with  electrons.  Appendix  M 
describes  the  geometry  models. 


The  energy  is  deposited  by  subroutine  ELOSS.  It  used 


the  range  energy  described  by  Fraley  e_t  al_.  (Ref  9:47S), 
namely 


(114) 


where 

s  is  distance  in  centimeters 
Po  =  0.213  grams  per  cubic  centimeter  for  solid  DT 

and 

3.5M*V  (H5) 


Ea  being  the  residual  energy  of  the  alpha  particle  in  Mev. 
The  first  term  represents  energy  deposited  in  electrons  and 
the  second  term  represents  energy  deposited  by  ions.  The 
fraction  of  energy  partitioned  to  the  ions  .  Fi.  during  alpha 
particle  transit  is  described  by 


D 


_  (3*  - 


a 


-i 


(116) 


The  energy  deposition  rate  is  much  higher  near  the  end 
of  the  alpha  particle  path  due  to  ionic  collisions  (Ref  9:475). 
Subroutine  ELOSS  models  this  by  comparing  the  electron  and 
ion  terms  and  stopping  the  particle  in  the  current  zone  when 
the  ion  term  dominates. 

The  range  of  the  alpha  particle  when  only  the  electron 
term  o'f  Eq  (114)  is  considered  is 

(y'jo.oa^  i  *  o.nAjj'  (0']J 

The  range  when  only  the  ion  term  is  considered  is 

JU  *  (j-) ‘5[ 1  * °-07J'in[T'’4(j)V1]j 

These  relationships  are  used  in  subroutine  FINDDS  dividing 
Xu  by  50  to  attain  an  integration  increment.  When  the 
electron  temperature  is  below  20  kev,  Xag  is  used  denoting 
electron  collision  domination.  Above  20  kev,  X  .  is  used 
to  model  ion  collision  domination. 

This  integration  increment  is  always  computed  prior  to 
energy  deposition.  The  actual  integration  is  accomplished 
in  ELOSS  using  a  4th  order  Runge-Kutta  method. 


I 

(117) 


(118) 
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Subroutine  HTFLX 


Subroutine  HTFLX  calculates  the  energy  transfer  between 
radiation,  electron  and  ion  components.  Flux  continuity 
at  cell  boundaries  is  used  in  electron  and  radiation  diffusion 
calculations. 

Constants  used  include  pi,  electron  mass  and  charge,  the 
speed  of  light,  Boltzmann's  constant,  the  fine  structure 
constant,  and  the  Stephan-Boltzmann  constant.  Constants  for 
coefficients  of  thermal  conductivity,  equilibration  frequency, 
Compton  coupling  and  bremsstrahlung  coupling  coefficients 
are  grouped  in  front.  Constant  terms  from  the  Coulomb  log¬ 
arithm  are  found  in  initial  subroutine  stages.  Heat  and 
radiation  fluxes  are  initialized  to  zero. 

Electron-ion  energy  exchange  is  the  first  process 
addressed.  The  Coulomb  logarithm  for  electron-ion  energy 
exchange  is 


£r\  _A*e  c  “  ^ 


(84) 


The  minimum  value  of  the  Coulomb  logarithm  is  set  at  1. 
The  term  A/(ZN  p)  is  coded  as  volume  per  number  of 

a 


electrons . 


Using  this,  the  thermal  conductivity  is  given  by  (Ref  26:144) 


»  2  i- 

Ke'40\Try  e4^/1 


(19) 


The  factor  e  is  due  to  the  secondary  electric  fields  origin 
ated  by  the  flow  of  charged  particles.  The  factor  6T  is 
dependent  on  the  Z  of  the  material  and  is  included  to  model 
a  real  gas  vice  the  Lorentz  gas  on  which  the  derivation  of 
Eq  (19)  is  based  (Ref  26:144).  The  product  efiT  is  modelled 


0.43  2 _ 

(3.44  +  2  +  o.ztJnL*]) 


(83) 


The 


equilibration  frequency, 


is 


3 


(82) 


The  quantity 


( 


av4 


(119) 


where  N  .  is  the  number  of  electrons  in  zone  j 
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The  updated  electron  temperature  comes  from  a  combination 
of  ion  electron  energy  transfer  and  electron  heat  conduction. 
This  is 
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At 


By 

2  k  hit 


(Tt-70 


(123) 


This  equation  is  discussed  in  Appendix  N.  The  ion  temper¬ 
ature  is  then 


dli 

<u 


-  *ef  (Ti-T*) 


(124) 


where  again  Z  is  net  charge. 

The  electrons  are  much  lighter  than  the  ions  and  so 
conduct  much  more  energy  through  diffusion  than  do  ions. 
Consequently,  ion  heat  conduction  is  not  considered. 

Radiation-electron  energy  transfer  is  treated  using 
a  Rosseland  mean  free  path 

from  Fraley  ejt  al . 
coefficient 

Kr  = 


/■  3T, 


Vx 


(90) 


(Ref  9:476).  A  radiation  diffusion 


/6<T  n  -r 
-3-  *p  ,r 


(91) 


is  then  evaluated.  Radiation  flux  is  then  calculated  by 

RF LUX  =  r*Kr 

dr 


(125) 


79 


r 

r 


where  the  spatial  derivative  of  the  radiation  temperature 
is  treated  similarly  to  the  spatial  derivative  of  the 
electron  temperature. 

A  bremsstrahlung  Gaunt  factor  is  computed  using 


(126) 


where  Y  =  Tr/Tg  •  The  Gaunt  factor  for  bremsstrahlung 
processes  is  discussed  in  Appendix  G.  This  Gaunt  factor 
is  used  in  computing  the  bremsstrahlung  coupling  coefficient 
which  is 


f  e4N« 


(86) 


The  Compton  coupling  coefficient  is  denoted 


(87) 


These  are  added  to  arrive  at  the  complete  coupling  coeffi¬ 
cient  or 


(85) 


These  coupling  coefficients  are  discussed  in  Chapter  II 
and  also  in  Appendices  F  and  I. 
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Specific  heats  are  required  and  are 


3  hj 
*  m 


and 


Cvr  =  iS'T'Tr 
CP 


The  radiation  temperature  can  then  be  updated  using 

hL(r..-n). 


pCyfV, 


This  equation  is  discussed  in  Appendix  © 
temperature  is  updated  using 

iS.  =-  (Te'Tc) 

dt 


The  electron 


Because  of  the  mass  difference  between  ions  and 
electrons,  radiation  interaction  with  electrons  is  highly 
preferential.  Radiation- ion  energy  exchange  is,  therefore, 
not  considered. 


Subroutine  OUTPUT 

This  subroutine  can  be  called  at  timely  locations  to 
provide  a  variety  of  information.  A  summary  of  the  data 
printed  on  hard  copy  is  given  in  Tables  S  and  6. 
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Table  5 


Cyclic  Data  Printed  by  Subroutine  OUTPUT 


Parameter 

Dimension 

Cycle  Number 

Pure  Integer  Number 

Total  Time 

Nanoseconds 

Time  Step 

Picoseconds 

Input  Energy 

Kilojoules 

Output  Energy 

Kiloj  oules 

Table  6 

Zone  Data  for  the  Given  Cycle  Number 

Printed 

by  Subroutine  OUTPUT 

Parameter 

Dimension 

Zone  Number 

Pure  Integer  Number 

Radius 

Millimeters 

Cell  Wall  Velocity 

Millimeters  per 
nanosecond 

Reactions 

Reactions  per  time  step 

Density 

Density  times  normal 

Electron  Temperature 

kev 

Ion  Temperature 

kev 

Radiation  Temperature 

kev 

Output 

Joules 

Heating 

Neutron  Fluence 

Joules 
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The  cyclic  information  is  printed  at  the  top  of  the  page 
followed  by  zone  information  in  columns.  Neutron  fluence 
is  a  parameter  of  interest  currently  not  computed  in 
program  MOXNEX.  Conversions  are  made  in  subroutine  OUTPUT 
to  attain  the  desired  parameter  dimensions. 


IV.  Validation  Arguments 

Introduction 

Code  validation  for  the  MOXNEX  program  is  a  goal  that 
must  be  satisfied  to  enable  its  use  as  a  vehicle  for  further 
study.  Because  the  code  is  not  yet  running,  comparison 
with  published  data  is  not  possible.  Steps  were  taken, 
however,  to  analyze  specific  analytic  cases  and  to  validate 
program  modules.  The  efforts  to  verify  the  MOXNEX  coding  and, 
in  particular,  the  HYDRO  and  HTFLX  subroutines  are  discussed 
in  the  following  paragraphs. 

Analytic  Cases 

The  radiation-electron  energy  transfer  coefficient 
Aer  useci  in  subroutine  HTFLX  is  discussed  in  Chapter  III. 
Chapter  III  also  includes  a  comparison  between  the  radiation- 
electron  energy  transfer  coefficients  used  for  this  study 
and  those  coefficients  used  in  Reference  9.  Chapter  III 
further  points  out  that  the  coefficient  used  by  Reference  9 
is  dimensionally  incorrect.  It  also  notes  the  MOXNEX 
bremsstrahlung  and  Compton  coupling  coefficients  differ 
from  those  used  in  Reference  9  by  a  factor  of  ptq  .  This 
means  the  radiation-electron  energy  transfer  coefficients 
used  in  the  two  studies  differ  by  a  factor  of  Tq  * 

2.8178  x  10  ^  centimeters  (Ref  3:501),  or  a  difference  of 
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almost  13  orders  of  magnitude.  A  difference  of  this  size 
in  the  radiation-electron  energy  transfer  coefficient  should 
make  a  substantial  difference  in  the  results  attained  from 
the  two  codes  and  is  further  discussed  in  Chapter  V. 

Critical  examination  of  the  processes  that  are  included  in 
the  radiation-electron  energy  transfer  coefficient,  then, 
is  central  to  this  study  in  order  to  confirm  the  value  of 
this  coefficient.  Appendices  include  derivations  of  the 
radiation-electron  energy  transfer  coefficient  components, 
the  bremsstrahlung  coupling  coefficient  and  the  Compton 
coupling  coefficient.  The  locations  in  appendices  of 
derivations  of  the  bremsstrahlung  and  Compton  coupling 
coefficients,  along  with  other  derivations  and  analytic 
verifications  used  directly  in  MOXNEX  coding,  are  summarized 
in  Table  7. 

In  many  cases  supporting  calculations  for  derivations 
are  presented  in  prior  appendices.  Cross  section  derivations 
for  short  and  long  range  Coulomb  collisions  are  derived  in 
Appendix  C.  Appendix  E  discusses  cross  sections  for 
bremsstrahlung  and  inverse  bremsstrahlung  collisions  which 
support  bremsstrahlung  coupling  coefficient  calculations. 
Bremsstrahlung  Gaunt  factor  modelling  is  addressed  in 
Appendix  G.  The  enhancement  factor  for  bosons  used  in  both 
bremsstrahlung  and  Compton  calculations  is  derived  in 
Appendix  L.  The  average  energy  exchange  during  a  Compton 
collision  is  discussed  in  Appendix  H.  The  geometrical  models 
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Table  7 


Analytic  Verifications  of  Specific  Terms  and  Equations 
Used  for  MOXNFX  Coding 


Appendix 

Derivation 

Subroutine  Page 

Where  Used 

F 

Bremsstrahlung  Process 

Coupling  Coefficient 

HTFLX 

F 

Bremsstrahlung  Gaunt  Factor 

HTFLX 

I 

Compton  Process  Coupling 
Coefficient 

HTFLX 

J 

Lagrangian  Energy  Equation 

HYDRO 

K 

Equivalent  Fermi  Electron 
Temperature 

HYDRO 

L 

Radiation  Temperature  Adiabatic 
Update  Equation 

HYDRO 

N 

HTFLX  Electron  Temperature 
Update  Equation 

HTFLX 

0 

HTFLX  Radiation  Temperature 
Update  Equation 

HTFLX 

used  for  subroutine  ALPHA1  are  presented  in  Appendix  M. 

These  derivations  and  mathematical  discussions  are  presented 
to  support  presentations  in  Chapters  II  and  III,  but  also 
to  provide  an  analytical  check  on  the  self  consistency  of 
the  coupled  rate  equations. 

An  error  was  discovered  in  the  original  HTFLX  differencing 
scheme  for  temperature  update  equations  that  required  cor¬ 
rection.  Temperature  update  equations  in  subroutine  HTFLX 
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use  the  first  derivative,  dT/dr  ,_where  dT  denotes 


either  the  electron  or  radiation  species  temper atiire ,  This 
temperature  gradient  is  approximated  by  the  quotient  of 
finite  differences  between  cells.  The  temperatures  used 
for  this  finite  difference  must  be  the  temperatures  at 
the  time  of  the  subroutine  call.  The  temperature  differ¬ 
ences  in  the  original  equations  used  already  updated  temper¬ 
atures  in  the  outer  3  cells  yielding  an  erroneous  gradient. 

A  new  differencing  scheme  is  now  employed  to  ensure  that 
the  gradient  is  computed  in  all  cells  using  the  temperature 
values  at  the  time  of  the  subroutine  call. 

Other  discrepancies  were  found  in  the  signs  of  the 
temperature  update  equations  in  subroutine  HTFLX.  These 
sign  errors  were  changed  ensuring  that  energy  flows  in 
the  correct  direction. 

Quasistatic  Equilibrium  Study  for  Subroutine  HYDRO 

Subroutine  HYDRO  updates  species  temperatures  based 
only  on  adiabatic  changes.  No  energy  transfer  mechanisms 
are  modelled  in  the  subroutine.  Consequently,  species 
temperatures  should  be  updated  only  by  cell  expansion  or 
contraction  independent  of  differences  in  temperatures 
between  species.  If  the  species  temperatures  are  nonzero 
and  approximately  the  same  between  cells  within  the  species 
inside  the  microsphere,  cell  expansion  should  be  most 
rapid  in  the  outer  cells  since  the  species  pressures  outside 
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of  the  last  cell  are  very  close  to  zero  (all  species 
temperatures  outside  of  the  last  cell  are  initialized  to 
1.  x  10  kev).  After  a  time  period  short  compared  to  the 
time  required  for  microsphere  kinetic  disassembly  and  in 
the  absence  of  other  processes,  the  inner  cells  should 
demonstrate  no  temperature  changes  or  energy  exchanges 
and  the  outer  shells  should  have  decreased  temperatures 
due  to  expansion. 

These  phenomena  were  tested  and  observed  by  calling 
subroutine  HYDRO  consecutively  100  times  after  initializing 
species  temperatures  at  equilibrium  conditions  and  then  at 
various  non-equilibrium  conditions.  Using  40  cells,  no 
temperature  changes  were  observed  in  the  inner  34  cells  in 
any  calculation  and  temperatures  lower  than  initial  temper¬ 
atures  were  seen  in  at  least  the  outer  3  cells  in  all 
species  in  all  calculations.  Table  8  lists  the  initial 
temperatures  and  final  temperatures  in  cell  40  and  the 
total  time  for  the  run  for  the  4  runs  completed.  The  reason 
that  total  time  on  the  second  run  is  different  from  the 
other  three  is  that  a  different  time  step  was  used. 
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Table  8 


Quasiequilibrium  Studies  for  Subroutine  HYDRO 
(100  Iterations,  All  Temperatures  in  key) 


Electron  Temp  Ion  Temp  Radiation  Temp 

Run  Initial  Final  Initial  Final  Initial  Final 

No.  (All  Cells)  (Cell  40)  (All  Cells)  (Cell  40)  (All  Cells)  (Cell  40)  TOTAL  TIME 


1 

6 

2.770 

6 

2.770 

6 

4.954 

0.0808 

nsec 

2 

6 

2.647 

3 

1.323 

6 

4.898 

0.1 

nsec 

3 

6 

2.796 

6 

2.796 

3 

2.483 

0.0808 

nsec 

4 

3 

1.505 

6 

3.010 

6 

5.056 

0.0808 

nsec 

Point  Explosion  Calculation  for  Subroutine  HYDRO 


An  analytic  solution,  due  to  Sedov  (Ref  25)  is  available 
for  pressure  at  a  shock  wave  front  resulting  from  an  explosion 
at  a  dimensionless  point  (Ref  25:238-240).  This  solution  for 
spherical  geometry  is 


p  -A  2y-  (v-fy 

5>'  % 


(131) 


where 

o  5  U/  A 


(132) 


where 
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(133) 


(134) 


In  these  equations 

PQ  is  ambient  pressure 

Ps  is  pressure  at  the  shock  front 

/  is  the  ratio  of  specific  heats 

rs  is  the  distance  of  the  shock  front  from  the 
sphere  center 

Ec  is  the  energy  of  the  explosion  charge 

The  constant  is  a  function  of  /  and  is  approximately 
0.84  for  ^  =  1.4  and  0.48  for  /  =  1.67  (Ref  25:231-240). 

A  hydrodynamics  scheme  may  be  validated  by  initiating 
an  explosion  obtaining  values  for  pressure  at  the  shock 
front  and  shock  front  distance  from  the' sphere  center,  and 
comparing  these  values  with  values  obtained  from  the  analytic 
solution.  A  dimensionless  point  can  only  be  approximated 
in  a  hydrodynamics  scheme  by  placing  the  charge  energy  in 


90 


the  smallest  possible  volume.  This  should  lead  to  error 
in  the  form  of  small  random  deviations  from  the  analytic 
solution. 

The  calculations  were  done  using  =1.67  and  compared 
to  the  Sedov  solution  curve  (Ref  27:142)  for  =  1.4  .  The 
slope  and  shpae  of  this  curve  was  duplicated  in  each  calcu¬ 
lation.  The  differences  in  the  ratios  of  specific  heats, 
i  ,  anc.  the  scaling  of  EQ  with  ^  is  felt  to  be  the 
reason  that  actual  coincidence  was  not  attained.  A  solution 
for  %  =  1.67  can  be  constructed  using  Eqs  (131)  through 
(134).  It  is  felt  that  coincidence  with  such  a  solution 
would  be  attained  with  a  point  explosion  calculation  in 
subroutine  HYDRO.  These  calculations  were  not  pursued, 
however,  in  view  of  limited  time  and  resources. 
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V.  Recommendations 


Code  Completion 

Without  code  completion,  the  work  presented  is  largely 
inconclusive.  Code  completion  and  subsequent  comparison 
with  published  data  is  central  to  any  further  code  develop¬ 
ment  or  modification. 

Time  Step  Selection 

Currently,  the  time  step  used  for  the  entire  MOXNEX 
program  is  determined  in  subroutine  HYDRO.  As  noted  in 
Chapter  III,  this  time  step  is  defined  by  a  fraction  (less 
than  or  equal  to  1/2)  of  the  time  required  for  a  shock  wave 
to  travel  across  the  thinnest  microsphere  zone.  Using  this 
time  step,  erratic  values  of  the  temperature  gradient  are 
observed  in  subroutine  HTFLX  resulting  in  temperature  changes 
of  as  much  as  three  orders  of  magnitude  in  electron  tempera¬ 
ture  during  one  iteration.  These  dramatic  variations  in 
the  temperature  gradient,  as  well  as  the  resulting  electron 
temperatures,  are  uncharacteristic  of  inertially  confined 
fusion  microspheres  (Refs  4,  6). 

Consideration  may  be  given  to  redefining  the  time  step 
selection  method  to  obtain  time  step  values  yielding  gradient 
values  compatible  with  HTFLX  equations.  One  method  that 
could  be  employed  would  be  modifying  the  program  to  allow 
the  electron  temperature  and/or  the  radiation  temperature 
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update  modules  to  define  the  time  step  for  the  entire  program. 
The  first  derivative,  dT/dr  ,  could  be  iterated  until  a 
favorable  value  of  dt  is  attained.  A  favorable  value  of 
dt  would  be  the  largest  value  of  dt  that  yielded  well 
behaved  temperature  changes  between  cells.  This  value  of 
dt  could  then  be  used  as  the  time  step  for  the  entire  program. 
This  process,  however,  may  yield  a  very  small  time  step.  This 
time  step  may  result  in  stable  difference  equations  in  sub¬ 
routine  HTFLX,  but  also  insignificant  changes  in  other  sub¬ 
routines.  For  instance,  cell  walls  in  subroutine  HYDRO  may 
move  so  little  that  their  movement  is  barely  consequential 
so  that  an  increased  number  of  iterations  is  required  for 
the  overall  run.  A  very  small  time  step,  then,  may  greatly 
increase  the  total  number  of  computer  operations  required  per 
run  and  inefficiently  treat  the  overall  problem. 

Efficiency  might  be  regained  by  letting  the  HTFLX  and 
HYDRO  subroutines  proceed  at  separate  time  steps  character¬ 
istic  of  the  processes  modelled.  If  the  characteristic  time 
step  for  the  HTFLX  subroutine  is  significantly  smaller  than 
the  time  step  computed  in  the  HYDRO  subroutine ,  HTFLX  could 
be  sub-cycled  independently  from  the  rest  of  the  program. 

When  total  time  over  many  HTFLX  iterations  is  equal  to  the 
HYDRO  time  step,  the  HTFLX  subroutine  can  be  interfaced  into 
the  rest  of  the  program.  This  approach  may  provide  stable 
equations  throughout  the  MOXNEX  code  and  still  achieve  a 
minimum  of  computer  operations. 
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Once  the  code  is  completed,  validation  by  comparison 
with  published  literature  should  follow.  Brueckner  and 
Jorna  give  graphical  data  of  ion  temperature  versus  micro¬ 
sphere  radius  at  different  times  for  40  micrometer  deuterium- 
tritium  microsphere  with  a  uniform  initial  density  of  600 
grams  per  cubic  centimeter  and  initial  uniform  electron  and 
ion  temperatures  of  5  kev  in  the  center  six  micrometers 
(Ref  4:361)  to  which  MOXNEX  results  may  be  compared.  This 
data  is  seen  in  Figure  15.  The  MOXNEX  model  incorporates 
many  features  of  the  Brueckner- Jorna  code.  The  thermo¬ 
nuclear  reaction  rate  modelling  for  MOXNEX  is  taken  directly 
from  the  Brueckner- Jorna  coding  description  (Ref  4:358). 

Some  differences,  however,  exist  in  the  equations  used 
for  charged  particle  and  neutron  heating,  radiative  transfer 
and  the  equation  of  state.  A  tabular  Fermi-Thomas  equation 
of  state  is  used  for  the  electron  species  (Ref  4:357)  which 
should  model  electron  energy  much  more  closely  than  does  the 
MOXNEX  coding  which  accounts  for  electron  degeneracy  only  at 
low  temperatures  and  uses  a  perfect  gas  law  as  the  equation 
of  state.  Still,  similarities  in  the  processes  modelled 
lead  to  anticipation  the  MOXNEX  code  will  duplicate  the 
curves  seen  on  page  361  of  Reference  4  by  much  less  than  an 
order  of  magnitude  for  the  same  input  parameters. 
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Figure  15.  The  propagation  of  the  burning  front  in 
initially  uniformly  compressed  DT.  The  initial 
density  is  600  gm/cm3  and  the  initial  radius 
40  pm.  The  central  region  initially  heated  to 
5  kev  has  a  6  pm  radius.  The  curves  are  labeled 
by  the  time  in  units  of  10" 10  sec  (Ref  4:361). 
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The  radiation  temperature  used  by  Brueckner  and  Jorna 
are  not  given  for  the  above  example.  The  radiation  tempera¬ 
ture  in  the  MOXNEX  program  might  be  initialized  at  1  kev 
following  the  lead  of  Reference  9,  with  the  hope  of  dupli¬ 
cating  the  Brueckner- Jorna  data. 

A  second  basis  paper,  written  by  Gary  S.  Fraley  and 
others,  gives  extensive  data  for  numerous  microsphere  sizes, 
densities,  and  temperatures  on  pages  479  through  483  of 
Reference  9.  The  MOXNEX  coding  of  alpha  particle  heating 
is  patterned  after  equations  appearing  in  this  paper.  The 
coding  of  the  energy  transfer  processes  in  MOXNEX  for  species 
temperature  update  in  the  three  temperature  model  closely 
parallel  the  equations  of  Appendix  B  of  Reference  9.  A 
tabular  Fermi-Thomas-Dirac  equation  of  state  is  used  in  the 
Fraley  code  (Ref  9:474),  so  once  again  differences  may  be 
anticipated  when  MOXNEX  results  are  compared  with  Fraley  code 
results  due  to  differences  in  the  equation  of  state. 

A  large  difference  exists,  however,  between  the  MOXNEX 
code  and  the  Fraley  code  in  the  values  of  the  radiation- 
electron  energy  transfer  coefficient.  This  difference, 
discussed  in  Chapters  III  and  IV,  may  result  in  a  much 
different  energy  partitioning  among  species  between  the 
two  codes.  This  energy  partitioning  difference  will  affect 
many  other  areas  which  are  temperature  dependent  such  as 
thermonuclear  reaction  rates  and  species  pressures.  The 
disagreement  should  propagate  and  magnify  as  the  codes 
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progress  in  time.  The  degree  of  similarity  between  the  MOXNEX 
code  and  the  Fraley  code  hinges  on  the  relative  importance 
of  the  radiation-electron  energy  transfer  process.  This 
process  may  be  dominated  numerically  by  thermonuclear,  Cou- 
lombic,  or  other  processes  in  both  codes  in  which  case  the 
MOXNEX  results  should  closely  approximate  the  Fraley  results 
because  of  the  similarities  in  the  processes  modelled  and 
similarities  in  the  equations.  Much  better  than  order  of 
magnitude  agreement  may  be  seen. 

If,  on  the  other  hand,  other  processes  do  not  numerically 
dominate  the  radiation-electron  energy  transfer  process, 
correlation  between  the  MOXNEX  code  and  the  Fraley  code 
should  be  hard  to  recognize.  Radiation-electron  energy 
transfer  should  be  a  significant  process  (Refs  4,  10,  18,  26). 
Doubt  exists  whether  this  process  is  so  totally  dominated  that 
a  coefficient  difference  of  almost  13  orders  of  magnitude 
will  be  inconsequential. 

Nevertheless,  Reference  9  reports  agreement  with  other 
published  literature  (Clarke,  J.S.,  H.N.  Fisher,  and  R.J. 

Mason.  "Laser-Driven  Implosion  of  Spherical  DT  Targets  to 
Thermonuclear  Burn  Conditions,"  Physical  Review  Letters,  30 
(2):  89-92,  (January  1973).)  in  total  microsphere  energy 
yield.  The  three  temperature  code  used  in  this  paper  is  the 
basis  for  the  three  temperature  model  used  for  the  Fraley 
code  (Ref  9:474).  Both  Reference  9  and  the  paper  by  Clarke 
and  others  were  written  at  the  Los  Alamos  Scientific  Laboratory 
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between  1972  and  mid-1973.  In  fact,  one  author,  R.J.  Mason, 
is  a  common  author  for  both  papers.  The  possibility  exists 
that  the  agreement  in  results  is  at  least  in  part  a  conse¬ 
quence  of  the  similarity  in  their  coding.  Further,  the 
possibility  exists  that  the  radiation-electron  energy  trans¬ 
fer  coefficient  used  in  both  codes  is  the  same. 

With  this  background  in  mind,  the  Brueckner- Jorna  code 
results  should  be  considered  as  the  primary  comparison 
standard.  Differences  in  results  between  the  MOXNEX  code 
and  the  Fraley  code  might  be  explained  by  the  difference  in 
the  radiation-electron  energy  transfer  coefficient.  Results 
from  both  codes  are  computed  using  one  dimensional,  spherical 
Lagrangian  coding. 

Another  inertial  confinement  fusion  code  to  which  MOXNEX 
results  could  be  compared  are  the  LASNEX  code  results  from 
the  Lawrence  Livermore  National  Laboratory.  LASNEX  models 
sophisticated  physics  in  a  large  and  expensive  computer 
simulation  to  support  experimental  laser  fusion  studies  at 
Lawrence  Livermore  National  Laboratory.  Some  LASNEX  code 
results  might  be  attained  from  recent  volumes  of  the  Lawrence 
Livermore  National  Laboratory  laser  program  annual  reports. 
Additionally,  these  annual  reports  contain  information  con¬ 
cerning  experimental  laser  fusion.  A  strong  validation 
argument  is  made  if  MOXNEX  results  agree  with  LASNEX  results 
or  experimental  results  for  total  yield,  disassembly  times. 
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neutron  fluence,  and  species  temperature  profiles  in  time 
within  an  order  of  magnitude. 

Driver  Incorporation 

A  final  recommendation  is  to  incorporate  an  energy- 
delivery  mechanism  or  driver  into  the  MOXNEX  code.  Driver 
incorporation  will  add  completeness  to  the  MOXNEX  code  and 
make  it  more  comparable  to  other  studies  of  its  type,  for 
instance,  the  Brueckner- Jorna  code.  Study  may  be  augmented 
in  many  additional  areas  such  as  laser-plasma  coupling, 
particle  beam-plasma  coupling,  and  microsphere  design  by 
using  and/or  further  developing  the  MOXNEX  code.  Procedures 
outlining  driver  incorporation  to  a  fusion  microsphere  are 
included  in  the  paper  by  Brueckner  and  Jorna  (Ref  4). 
Additional  modelling  information  can  be  found  in  recent 
laser  program  annual  reports  by  the  Lawrence  Livermore 
National  Laboratory. 
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Program 


Listing 
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i  l  •">= 

1 30= 
14  0= 

01> 


rt-  M0>'NE>' 

r*  I  ME  '5  T  ~ri  R  1 0 1  >  .  AMASS  (101).  VEL  1  ■ 

iMHEf-  SION  H  '  1  0 1  •  -HR'  101  .«  .  PUSH  <  101  >  .  ELEC  U  01  >  -TEClOl )  . 

*  Till 0 1 >  ,  TR  vK.;)  .  1  - 1  tC  i  '  .  - 1"  ■'  1  r>  1  n  '•■'Ail 

DIMENSION  RELRHO ( 101 ). THEAT (101 >. DELV •  • 


I SO=Ct**Y  ******** * * ***  * *** * Y  ******** ******* 

160=C  HYDRODYNAMICS  VARIABLES - 

1‘  0=C  R(l)  IS  RADIUS  DR  ITH  ZONE  IN  CM 

1S0«C  AMASS (I)  IE  TOTAL  MASS  OR  ITH  ZONE  IN  ERA!  S 

l?OC  F ( I )  IS  TOTAL  PRESSURE  IN  ITH  ZONE  IN  EPS  "CM3 

2' >o=C  0(1)  IS  ARTIFICIAL  VI3CC3ITV  IN  ITH  ZONE  IN  ERG/CM3 

2 1 0=C  VEL  < 1 >  IS  VELOC  IT  Y  OR  ITH  ZONE  IN  CM /SEC 

220=C  V(I)  IS  VOLUME  OF  ITH  ZONE  IN  CM3 

23  =C  A'  *  Y  A  *  Y  w  A  *  *  *  *  #  Y  *  ***  *  *  A  Y  *  Y  A  *  Y  k  *  A  A  *  *  *  *  *  A  Y 


24o=C  INDIVIDUAL  PARTICLE  VARIABLES - 

250=C  H ( I >  IS  TOTAL  NUMBER  C"  D+T  ATOMS  IS  ITH  ZONE  (50-50  MIX) 

260=C  HE (I)  IS  TOTAL  NUMBER  OF  ALPHAS  IN  ITH  ZONE 

270=C  PUSH (I)  IS  TOTAL  NUMBER  CF  PUSHER  "ATOMS"  IN  ITH  ZONE 

2S0=C  ELEC  (I)  IS  TOTAL  NUMBER  OF  ELECTRONS  IN  ITH  ZONE  (NEUTRAl.) 

270 =  C  TE  ( I )  IS  ELECTRON  TEMP  I!  JTH  ZONE  IN  KEV 

300= C  T I  ( I >  IS  ION  TEMP  IN  ITH  ZONE  IN  KEV 

310=C  TR ( I >  IS  RADIATION  TEMP  IN  ITH  ZONE  IN  KEV 

320=C  FEAC(I >  IS  TOTAL  NUMBEF  OF  REACTIONS  IN  ITH  ZONE  THIS  LT 

330=  C***  *  Y.  Y  Y  *  %%.%%  *  %  *  Y  Y  i  *  :i  t  Y  Y  Y  t.  Y  Y  Y  Y  Y  *  *  t  Y  Y  Y  Y  v 

340=C  SINGLE  VARIABLES - 

350=C  DT  IS  TIME  STEP  IN  SECONDS 

360---C  NSPACE  IS  NUMBER  OF  I  Or  ICS  IN  ENTIRE  PROBLEM 

370=C  NFUEL  13  NUMBER  OF  LAST  FUEL  ZONE 

380-C  ZPUSH  IS  AVERAGE  ATOMIC  '  •  BEK  CF  PUSHER  "ATOMS" 

390=0  Y*******  *  Y  *  Y  Y  Y.  Y  *  Y  Y  t  Y  Y  Y  Y  ********  Y  Y  Y  Y  Y  Y  Y 

4  00=C  T  ASKS - 

410=0  SET  INITIAL  CONDITIONS 

420=  CALL  GDATA  (F;,  DT.  TE.  TI .  TR,  NSPACE,  NFUEL,  TIME.  H,  HE.  PUSH 

. CYCLE. 

4  30=  1  ZF  UGH ..  PM  ASS ,  AMASS .  TMABS .  AMASS ,  ELEC .  V .  HEMASS .  DELV  ,FR.  HD .  P 

USHD. 

440=  Y  HMASS.  DFVENG.  TREAC.  TKF.AT,  VEL,  ALPHA) 

450=  DF;ENG  =  0.0 

460=  DO  10  1=1. NSPACE 

470=  READ (I)  =  0.0 

460=  RELRHO ( I )  =  1. 

490=  10  FLU ( I )  =  0.0 

500=  CALL  OUTPUT (R. CYCLE. TIME.  DT, DRENG.  NFUEL. NSPACE. VEL. RHO 

Z.RMASS 

510=  */..  ZF'USH.  ELEC.  TE,  TI .  TR,  REAC.  TREAC,  TREAT,  FLU,  FELRHO .  H,  HE.  PUSH. 

DRVENG ) 

520=C  UPDATE  HYDRO  VARIABLES.  PARTICLE  VARIABLES  BY  HYDRO  (WADE/W 

EBER) 

530=  CALL  HYDRO  (R. VEL . DT , TR. TI , TE. NSPACE . TIME. CYCLE, I BURN, 

PELF.  HO. 

540=  1  H.  HF  .  ELEC .  PUSH .  AMASS.  PMASS .  DMAS3,  TMASS .  HEMASS ,  PR .  DELV .  HD .  P 

USHD. 

550=  2  NFUEL. V) 
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CALL 


Z  ,  PM/  S3 

C.-7, 


UT  -ft.t'  CLi.T  Ji 


i  rUE: ..  N3PY .f„E  .  -  Ei  .  RHO 


ZPUSH,  ELEC.  TE.  TI  .  TF;  F U ' r  TREAT  TFrA~  H  ■’  r'1_  P'-»n.  U  ->E  P'  EH. 


5S0=C 

590- 


600= 
Z. PMASS 
610= 
DRVENG) 


s>  4 0= 

7  PMASS 
650= 

D  WENS ) 
660=C 
670= 
630= 
690= 
700= 

Z . PM ASS 
7 1 0= 
DRVENG ) 
720=C- 
730= 
MASS. 
740= 


CONFUTE  THERMOMUC  LEAF  REACTION  RATES,  REMOVE  DT  MASS 

CALL  TNBURM (NSPACE.  NFUEL.  DT.  V.  H,  T I . REAC . TREAC , I  BURN) 

CALL  OUTPUT (R. CYCLE  .  TIME. DT. DRENG. NFUEL. NSPACE. VEL, RHG 

ZPUSH, ELEC.  TE .  T I  «  TR.  REAC .  1  REAC .  1  HE  A  .  PLU .  F:ELF HO .  H .  HE, FLSH, 

DEPOSIT  NEUTRON  ENERG'- .  COUP  JTE  -1EUTRQN  OUTPUT 

Cf 'LL  NEUHEt-  i  (NS*-f  r  t:  .  k«  F>EA-._,  AL  FFh.  THEa  '  .  -i,  Hc^"> 

CALL  OUTPUT  (R.  CYCLE.  TIMS  .  DT.  DRENG.  NFUEL  .  Nr  3 ACE.  VEL. PH C 

7. .  CRUSH .  ELEC , TE , T 1 . TR .  REAC ,  TREAT  .  THE  AT  .  FLU ,  RELF  HD .  H .  HE .  PUS- . 

DEPOSIT  ALPHA  ENERGY.  ADD  ALPHA  MASS 
IALPH  =  0 

CALL  ALPHA 1 (R, AMASS . V . H .  HE .  F USH , ELEC , TE . T 1 . REAC . DT , NSRACE . 

+  NCULL.  2 PUSH .  PMACS  .AL PH.  TREAT..  ALPHA  '• 

CALL  OUTPUT  (ft.  CYCLE.  T  j  ME.  D!  ,  L  REMS.  NFUEL.  NS  PACE,  ■  Et. .  R--C 

*/.,  ZFUSH.  ELEC.TE.TI .  TR.  REAC,  TREAC.  THEAT.  FLU.  RELRHO.  H,  HE,  PUSH, 

CALCULATE  RADIATION  ENERGY.  HEAT  CONDUCTION 

CALL  HTFLX (R. TE,  TI ,  TR,  ELEC,  V,  H,  HE. PUSH. AMASS , NSRACE . H 

*  HEMASS.PMASS.DT.FR) 

CALL  OUTPUT < R . CYCLE .  T I ME .  DT  .  DRENG .  NFUEL . NSRACE . VEL . F HO 


760=  7.,  ZPUSH,  ELEC.  TE,  T  1 ,7ft.  REA(  ,  TREAC,  THEAT,  FLU,  RELRHO.  H.  HE,  FUEH. 

DRVENG) 

770=C'i  *4  ##*4**tt*#*  *4  *.**#**.♦**  t  ****.&**  *  tttt 

760=  END 

*’SO=C  ttt  tt*tt  if  ttt  t '  *  4  4  »  i*  4  *  i  4  4  .*  4  4  .f.  4  4*  •  *  *  4 

500=  SUBROUTINE  GDATA  (R. DT. TE.  1  1  .TP. M SPACE. NFUEL, TIME. H. HE. PUSH 

.CYCLE. 

.  .  L 1ST  300.  250)0 


000= 
. CYCLE, 
810= 
U5HD. 
620= 
B30= 
H  ( 1 0 1 )  , 
340= 
850= 
860'= 
870= 


SUBROUT IMF  GDATA  (P. DT.  TE, TI  .  TR.  NSPACE.  NFUEL, TIME. H, HE. PUSH 

1  ZPUSH . PMASS . DMASS , TMASS . AMASS , ELEC . V , HEMASE , DEL  V . PR , HD , F 

4.  HMA3S .  DRVEi-JS .  TREAC .  “'  MEAT .  VEL ,  ALPHA ) 

DIMENSION  FYlOl) . TE  < 101 > .TI (101) .TR(101) ,H(101 ■  . HE ( 101 ) .PUS 


340=  1  ELEC  (101).  AMASS  ( 1 0 1  )  .  DELV  ( 1  0 .!. )  .  FT-  (101)  ,  V  ( 1  0 1 )  .  VEL  ( 1 0 1 ) 

850=  DIMENSION  DRVENG (101 )  .  TREAC  < i 0 1 ) .THEAT- 101) . ALPHA (101) 

860'=-  INTEGER  CYCLE 

370=  PIE  =  4 . 4ATAN (  1 .  ) 

380=C  SET  MASSES  OF  COMPONENTS  <GMS) 

SF 0=  DMASS  =  3.34E-24  . 

900=  TMASS  =  5.0)27-24 

910=  HMAES  =  (PMASS  +  "MASS)  2. 

920'=  HEMA55=6. 64E-24 


960  s* 

970=C 

«S0= 

C-C 

i'  '00 =C 
1  ic- 


jl  i :  -■ 

: or o-c 
ios  •= 

1070-C 
1 05 0— ■ 
1090= 

1 100=C 
'  i  i  r.= 

i  1 20= 


••  "  ■  j_.  :  .  fcc.-*  t 

JMT-T  '  ~  '  'TV  P  >  T  7’sr 

TIME  -  0.0 

CYCLE  =  0 

BET  1  UMBER  OF  CELLS 
N5F  TCE=40 

NSF-U-NS-'V'Ct  M 

SET  HYDROSEN  DENSITY  (GM/CM3) 

HD~  ‘-■00 . 

TFT  ~  ~  jr  UIDTH  C1‘  AND  3  T  EM---'  . 

I  -:-  .  0001 
R'^7 

MINIM-  N  TIME  STEF 
DT= 1 . E- 1 2 

FIND  AVERAGE  Z  OF  PUSHER  ATOMS 
Z PUSH-- SORT  (37.  ) 

DO  1  N=i , NSFP1 

COMPUTE  WALL  POSITIONS  (CM)  AND  VOLUMES  (CM3) 
R(N)  =  DP  M t N— l > 

RP3  =  <F>  (N)+OR)  tt3 


1  130= 

V  ( N ) 

=  4.  ♦PTE *  ( RP3-PM3 )  / 

i  L  40= 

RM3  = 

3 

t.n 

ii 

n 

COMPUTE 

MAS 3  IN  CELL  (GMJ 

1 1  b  0  - 

hMASS 

( N )  =  HD¥V (N) 

i  170=0 

COMPUTE 

NUMBEF  OF  PARTICLES 

1 130= 

H  (N) 

=  AMASS (N) /HMA5S 

1 1 90= 

HE  <  N ) 

-  0 . 0 

1200=  PUSH  <M>  =  0.0 

1210=  ELEC (N)  =  H (N) 

1220=  ALPHA (N)  =  0.0 

1230=C  SET  VELOCITIES  (CM/SEC)  AND  VOLUME  CHANGES  (CM3) 

1240=  VEL(N)  =  0.0 

1250=  DELV (N) =0 

1260=0  INITIALIZE  RADIATION  PRESSURE 
1270=  PR  <N> -0. 

1280=0  SET  CELL  OUTPUT  (ERGS).  REACTIONS.  AND  CELL  HEAT  INS  (ERGS 
. 290=  DRVENG ( N ) =0 . 

1 300=  TREAC ( N ) =0 , 

1310=  THEAT (M>  =0. 

1320=C  AMBIENT  TEMPERATURE  FOP  CELLS  WITH  NO  ENERGY  INPUT  (KEV) 

1330=  TE ( N )  =  l.E-7 

1340=  TKN)  =  l.E-7 

1350=  TR(N)  =  l.E-7 

1360=  1  THEAT (N)  =  0.0 

370=C  YOU  HAVE  TO  HAVE  MORE  P-VALUES  THAN  OTHER  VAF I ABLER, 
13S0=C  SO  THAT  THERE  IS  A  VALUE  OF  RUH-1)  FOR  THE  LAST  ZONE 
1390=C 

1 400=0  SET  NUMBER  0^  FUEL  CELLS 
1410=  NFUEL  =  NS^ACE 

1420=C  INPUT  ENERGY  (KEV) 

1430=  CO  2  N=1 , NSPACE 

1440=  TR(N)=1. 

1450=  TKN)  =5. 

146  0=  2  TE  <N) =5. 

1470=  RETURN 


-L  T  '  .  r-  E1’  w  .  1  RE  AC  .  1 SUB  ■ . 
.  .FErCMOi  ’ 


A:.!.  TEi-L-  in  this  pt  init 


!  I  UN  I  ir  _  Rn  I  L’r 
-ITU-’  WITH  Jr  Jr  '  If 

INITIATED 
DO  2  I  "1 .  NFLIEL 

2  IF  ITI II ) . GE. TSTART)  IBURN=U 
I F  ( I BURN. t l 

CONFUTE  NUMBER  OF  REACTIONS  (D+T  ONLY)  IN  EACH  CELL  OF  FUEL. 
DO  3  1=1. MFUEL 
IF  ( T I ( 1 > . LT . 0 . 0 1 )  THEN 

reach:  i)  =  o.o 

GO  TO  3 

~r  T'.rr 


SIGVDT=  (3.41E-J4)  VTKO'JUEXR  (  <3.  o38*  TN13)  -  ‘  27 . 2  1 74  TNZ3  >  ) 

GO  TQ  5 

4  S J GVDT  =  (3.8E-i2)*TM2r  *EXR  <- 19. 02*TM13) 

5  FEAC  ( I )  =  DT*  (H  J>  ‘12’  xfilfiVBT/ (4. 0*V<  I)  ) 

REMOVE  BURNED  FUEL  N  >CLEI  FROM  CELL 

El '  I >  --  H  <  I  >  —  <  2  .  0*REAC  .[  '<  ' 

ADD  NUMBER  OF  REACT  JONS  IN  THIS  DT  TO  TOTAL  FOR  EACH  CEL.L 
TT.  CAT  1 1  '•  TREAC  •;  n  -  FEAC  <  I  > 

CONTINUE 

ADJUSTMENT  OF  ELECTRON  NUMBER  IN  EACH  CELL  WILL  BE  PERFORMED 

ANOTHER  SUBROUTINE  WHEN  THE  CHARGED  PARTICLES  PRODUCED  IN  T 

REACTION  HAVE  BEEN  TRANSPORTED  TO  OTHER  CELLS. 

UNITS  IN  THIS  SUBROUTINE  ARE"  TELE  IN  KElL  VOLUME  IN  CM  CUBED 

REACTIONS  IN  NUMBER  OF  REACTIONS  IN  EACH  CEL 

AND  SIGVDT  IN  PER  CM  CUBED  FEE  SECOND. 

RETURN 
END 

SUB. ROUT  1  NE  HYDRO  <  R. '  JtL .  DT  .TR.T1.TE,  N5RACE .  TIME ,  CYCLE  I  BURN . 

1  H,  HE.  ELEC.  FUSH.  AMASS.  F'MASS.  DMASS.  TMA3S,  HEMASS.  PR.  DELV.  HD. 

2  NFi'EL.V) 

D I  MENS I ON  R ( 1 OJ )  .  FT ( ] 01 > . PI ( 1 01 >  . PF (101), PR (101). 


■'  9'4« 

!  -f-O- 

1 980= 
1990= 
2000= 

203  0= 
2  4  0 — 

.  :>5<:  = 


0  ( lOl  .  OZL  iOi> .7F  (201  )  .Tt(l<  t  /  .  AH  C  5  •'  i  *1)  T1  .  !fl  •  . 
H  ( 101  )  „  HE  viol'.  P  !EH  (1>M,  .  E!  EC  v  1  .  1 2  . 

DEL11  ( 1 0 1  )  «  RELRHO  ( l  O  i  /  ,  RHO  (101/  .  ACC  <  1 0 1 ) 

T  .-A  PEEP  CYC_E 

SIGMA  =  5.67E-5  *  ( 1 . 1 o 1 L7 ) » *4 

A  =  4.  *  SIGMA  /  2.997E10 
HMAS3  =  ( DMASS  +  TMASS1/2. 

DTMAX  =  i.E-10 

CYCLE  =  CYCLE  +  1 
T IE  =  4 .  * AT AN  1 1 . / 

Du  3  0  I=L  .1  .SPACE-*- 1 


:F  •,  ]  BURN. ME. 


'B  (  ’  )  "H  ( I 


IS-rptJ  ~H  •:  ) 


•:<  ’SO  -  F  1  ( I  >  =  (H  ( I  >  ■‘•He. \  I  3  ■••PUSH  ( I )  /  tT I  <  i  )  ♦  1 . t-E-S  /'-■  ( I ) 

.  070=  F P  <  J  )  =  A  »  T*  (I) **4  /  3.. 

208 0=C  NOW  OBTAIN  EFFECTIVE  ELECTRON  TEMP  TO  ACCOUNT  “Of;  DEGENE 

2090=  EFFTt  =  0.00565*  (AMASS  (T  ■  / V  ( I  )  >  **  (2.  /3.  ) 

2100=  IF ( (I.GT.NFUEL) .OF . (EFFTE.LT.  TE ( 3 > > >  EFFTE=TE ( I > 

21 1 0=  PE ( I >  =  ELEC ( I )  *  EFFTE  *  1 . 6E-9  /V  ( I > 

0 120=  PT  ’ I )  =  PE ( I >  +  PI (I>  *  PR  ( I ) 

2 1 30= C  INITIALIZE  DENSITIES  AND  SET  0=0  EOF  FIRST  CYCLE 
2140=  IF  (CYCLE. EG. J  >  RHO < I > =AHA5F < I )  , 

.  i  5  •=  I E  ( CYCLE .  Eu.  X  )  0(1)  =0 . 

2160=  10  CONTINUE 

2170=  NOW  FIND  TIME  STEP 

2180=  DO  15  1=1. NSPACE 

2 1 90=  C  =  SORT  ( 5 .  /  3 .  #•  FT  ( I )  *.V  ( I  >  /  AMASS  <I>> 

2200=  FF. ACT  = .  5 

2  7  :  j  DT I  FRhCT $  ( E  (  i  *•  1 )  — F.  ( I  >  )  /C 

2220=  IF  (DTI .LT. DTMAX)  DTHAX=DT I 

2230=  IF  (DTMAX. LT.DT)  DT=DTMAX 

224  0=  1.5  CONTINUE 

2 ~ 5 0=C  NEXT  FIND  ACCELERATIONS  DC  CELL  WALL 

2260=  DO  17  1=1. NSPACE 

1270=  AN=- <PT ( 1  +  1 ) +Q( 1  +  1 ) — PT  < I ) —Q (I) ) 

2280=  IE  (I .  (ME  NSPACE )  AD=  <*(i+z>  -R  <r  > ) /l . *  ( PHD  <  1*1  +PHC  (I)  ) 

2290=  IF  (I.  ED.  NSPACE)  AD=  (F.  <  T-i  )  -R  ( I  )  *  (RHO  <1  +  1/  =RHO  ( I )  )  /2. 

2300=  ACC ( I ) “AN /AD 

231 G=C  NOW  FIND  VELOCITIES  OF  CELL  WALL 
2320=  VEL  ( I )  =VEL.  ( I )  +ACC  <  I )  *DT 

2330=  17  CONTINUE 

2340=C  FIND  NEW  POSITIONS  CF  CELL  MALI 
2350=  DO  19  1=1. NSPACE 

2360=  R  ( I  + 1 )  =F:  <  I  + 1 )  +VEL  ( I )  *  DT 

2370=  19  CONTINUE 

2380=C  CONFUTE  NEW  DENSITIES  AND  VISCOSITIES 
2390=  DO  20  1=1. NSPACE 

2400-  RHOOLD=FHO ( I ) 

2410=  VOLD  =  V(I) 

1420=  V(I)  =  4./3.*FTE*(R(I+l>**3-R<I>**3) 

2430=  DELV(I)  =  V(I)  -  VOLD 

2440=  RHO  ( I ) --AMASS  ( I ) /V  ( I ) 

24  5  0=  R  ELRHO  ( I  >  =F;HO  (I)/. 213 

24  60=  D I V=  .. 1 . -RELRHO ( I ) ) /DT 

14  0=  0(1 )  =RHO  ( I )  *  <R(  1+  1  >-F:  ( I )  /  *  'R  ( I  +  t)  -R  ( I  >  >  .fDIVtDIV 


108 


.  •  .  W1-.  .  UU>  U"  i  ) 

500=  Till)  =TI  I )  —  (2  .  /  i .  i  *  <F  1  ( I )  +u  <  1 )  (  i  ELV  <  i  }  /  •.  J  .  <  S— 9*  ‘.t-i  j  >  —  -.£  I 

)  )  ) 

..LIST  2500,2700600 


>  )  'i 

2510— 

..5*:  = 

2240= 
235.:  = 
2560=  30 
2570= 
2590= 

_ 5  v  o = 

2420= 

2470= 

2  44'I)= 

0  i  ) 

2450=C 
24  60=C 
247  >-r, 
24S'  »=C 
24 90 -C 
2oOO=C 

”  =  1  jsrr 

2  22  *'=C 
2530=C 


~T1 (1 /-  2 .  /  Z .  )  i  1  F  l  ( I  •  -t  0  k  I  )  )  *  DEL'v  ,  I  >  /  1  .  ©E-9  v  ( H  (  I )  +p» 

TE  (  t  •  =TF  M  )  -  .  2.  '7.  )  * -  E  (  7  )  tPELv  ■  I  ■  6  E-9  ELEC  <  I  )  > 

IP  1)  =  7 7-.'  ( I  *  •  ( vQl.D-BEL'.  I  7.  -'V  I , •  H . 25 

TF  >'  '  ;  7  >  ,  L7  .  Z.S--7)  71  I)  -1.F.-7 

“F-’Tr  I'  .  LT  i  E-~  T  :  =1  2-7 

IE< 3?  <  I  ,•  .LT.  l.C-7  tr(i  =l.E-7 

CONTINUE 

TIME=T2ME+DT 

RETURN 

END 

SUBROUTINE  1 1EUHEAT (MSP ACE.  P. REAC. ALPHr. TH^A^. V. H.HD 
P  T  HENS 7  ON  AL  PH A < 3 0 1 > . READ < 1 0 1 ) . EX ( 0 : 1 0 i '  . D  <  3 0 1 > . P ( 101) 

D I HENS I ON  THEA  T' 101) A  1 01 >  ,  H < 1 o) ) . HE ( 1 0 1 ) . AND (301'. AVGEI 


D,  T.  AMD  ME  CROSS-SECT  I  Oh'S  APE  f'-rpRQViHATF.LY 
.8  BARNS  FOR  14.1  MEV  NEUTRONS.  AN  AVERAGE 
ENERGY  PER  COLLISION  IP  FOUND  BY  AVERAGING  EACH 
REACTION’S  ENERGY  DEPOSITION  AND  THEN  SUMMING 
BASED  ON  NUMBER  FRACTIONS.  TREAT  RECEIVES  THE 
NENTRON  ENERGY  £4. F'G5 i  ~€t •  IK  EACH  CELL.  ALPHA 
IS  1  HE  EQUIVALENT  NUMBER  OF  3.5  MEV  ALPHA 
PARTICLES  FOR  THE  NEUTRON  ENERGY  DEPOSITED. 


254 0=C 

2550= 

SIGMA  =  . PE-24 

2560= 

TOT I NT  =  0. 

2570= 

DO  1  1  =  1,  NSF'ACE 

2530= 

BNO  =  H ( I )  +  HE ( I ) 

2590= 

AVGEN(I)  =  (  < .  5  *  H  ( 1 )  *  2 

35) /BNO) 


<  (.  5*H  ( 1  >  *  3  . 7625)  /  BNO) 


2600= 

2610= 

2620= 

2630= 

2640= 

2650= 

2660= 

2670= 

2680= 

2690= 

2700= 

2710= 

2720= 

2730= 

2740= 

2750= 

2760= 

FACE. 


*  ( (HE (1) *1.4) /BNO) 

ANO(I)  =(H(I)  +  HE (I))  /  V(I) 

TOT I NT  =  TOT I NT  +  REAC  ( I ) 

CONTINUE 
S  =  TOT I NT 
EX  < 0 )  =  1. 

DO  2  1  =  1 .  NSF'ACE 

F'.I  =  R(I  +  1)  -  R(I) 

EE  =  EXP  <— AMO ( I ) *SIGMA*RI ) 

EX  <  I )  =EE  *  EX(I-l) 

ACTIONS  =  S  *  (1.  -  EE)  *  EX (1-3) 

THEAT(I)  =  ACTIONS* AVGEN ( I ) * ( 1 . 6E-6 / 

ALPHA ( I )  =  THEAT(I)  /  <3. 5* < 1 ; 6E-6) > 

CONTINUE 

RETURN 

END 

SUBROUT  I NE  ALPHA  1 '( R .  AMASS .  V .  H .  HE  .  PUSH  . ELEC , T  E . 7  I ,  PEAC  DT  ,  NS 
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~C 

2800= 

01 ) , 
2810= 
2820= 


+  !  JR  1J '  _ .  _  '*  ;H  ,  r  M>  .  3  •,  r  H .  T1— r.A  i  AL  t-'ri-A 

D  T  MENS  I  ON  A!  -'HA  <  1  0 1  - 

D  3 HENS  3  ON  THETA  31)  F'LECd-' 1  £  .  u!  •!  101  > 

DIMENSION  R ( 101 ) , AMASS (101)  .‘O'  103  )  ,  HEdOl  >  ,  RUSH (101 ) . ELEC  ( 1 

+  TE  (101)  .  TI  (101  )  .REACdOl  >  .H(I03  /  .  THEAT  (101) 

DATA  THETA/O. , .7467380321.  .9689779077.  3  .  196752093.  1.3871656 


—  l)  ‘  IT 

:'950= 

’37*'.'= 


r  3 . 570796327.  3 . 75442~<  33.3. 94484056 

*•  2.  .  °46o462 1.3...  14  1592654/ 

DATA  PI /3. 3 41592654/ 

DO  10  1=1,1 0 1 
FE.LEC  ( I  >  =0  ■ 


153  61  4  746. 


2880= 

10  El ON ( I ) =0 .  ' 

2690= 

DO  2000  I 20NE=1  NFUEL 

2900= 

IF ( IALFH. EQ. 0)  GO  TO  250 

2910= 

DO  1000  I THET =1,11 

2920= 

ISTOF-O 

THET=THETA ( I THET ) 

2940= 

U=  1 . 

2950= 

X0= ( R ( I  ZONE ) +R ( I ZONE+ 1 ) ) /2 

2960= 

THETO=THET 

2970= 

INOW=I ZONE 

:980=C 

?990 =C*******  STAGE  FOR  BIRTH  SHELL 


3  I'QO-L 

3010=  IF ( ITHET. EQ.  1 ).  CALL  XYT6 ( INOW.  R,  XN,  THETN.  DI5T. I DUH . XO) 

3020=  IFdTHET.NE.  1. AND. THET. LT. PI/2.  >  CALL  XYT7  (XO,  INOW,  THETO,  R. 

DIST. 

"030=  I  CELL,  XN.  THFTJ) 

3040=  IF ( ITHET. EQ. 6)  CALL  XYT3 ( XO, XN. DIST . THETO. THETN. INOW, R> 

3050=  IFdTHET.NE.  1  1 -AND. THET. GT.F1/2.  )  CALL  XYT4  (XO,  XN.  DIET.  THEt 

0. 

3060=  +  THETN, INOW. R) 

3070=  IFdTHET.EQ.il)  CALL  XYT5  (XO,  XN,  DIST,  THETN.  INOW,  R> 

3080=  XQ=XN 

3090=  THETO=THETN 

31 00=C 

311 0=C COMPUTE  THE  INTEGRATION  INCREMENT 
ALF  HS=REAC ( I  ZONE ) / 1 1 . 

IF  (ITHET.  EQ.  1 )  ALF'HS=F:EAC  <  I  ZONE)  /  I  1 .  +ALPHA  (I  ZONE) 

TEMF'E=TE  <  I  ZONE) 

DENSE=AMASS ( I  ZONE ) /V ( I  ZONE) 

CALL  F I  HDDS  <  TEMF'E ,  DENSE .  DS .  DIST,  IDS) 

ICELL=I ZONE 

CALL  ELOSS  (DENSE.  TEMF’E,  U,  DS ,  D I  ST  .  EELEC .  E I  ON .  I  CELL.  I  STOP.  I  AL 

+  IDS, THEAT, ALPHS) 

IFdSTOF’.EQ.  3  )  ICELS=ICELL 
IF ( I STOP. EQ. 1 )  GO  TO  600 

3240=C******  STAGE  FOR  SUBSEQUENT  MOTION  OUT  OF  BIRTH  SHELL 
3250= C 

3260=  DO  500  I LOOF- 1.1000 


3 1 20= C 
3 1 30= 
3140= 
3150= 

3 1 60= 
3170= 
3190= 
3190= 
PH, 

3200= 

3210= 

3220= 

3230=C 
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IF  (  y0.  GF.  F  fNSPACF-M  )  '  GO  TO  600 

!  — *  Ci.  Tu  .  *.  -V.  'A "'6  i  T !Jii) >,P.  Ti^ETw.  DIET «  1CELL.  a‘.3 J 

I F  <  THE  TD ..  GT .  0 . 0 .  AND .  THETO  .LE.FI/2.  )  CAt_L  X  V  77  (  XO.  I  NON .  THE  ■  O 


+  R.D15T. ICELL. XN, THETN) 

I F ( THETO . GT . F I /2 . 0 . AND . THETO . LT . F I )  CALL  XYT8 ( I NOW . THETO . XO 


734  0= 
7300= 


33  30= 
PH. 

3390= 
3400= 
3410= 
3420= 
3430= 
3440= 
3450= 
3460= 
3470= 
3480= 
3490= 
3500= 
7510= 
3520  = 
3530'= 
3540= 
3550= 
3560= 
3570= 
3580= 
3590= 


THGOn.  XN.  DIET.  F:) 

IF  <  THETO.  EO.  PI )  CALL  XYT9  ( I  NOW.  XO  THETO.  F;.  ICELL,  DIST.  XN. 

IF  <  TCELL.  FT.  NSF'ACE)  GO  TO  3  000 
TEMPE=TE<  I  CELL  1 
DENSE=AMA9S ( T  CELL ) / V  < I CELL ) 

CALL  FINDDS  •  TEMF'E.  DENSE.  DS,  DIST,  IDS) 

CALL  ELOSS (DENSE. TEMPE. U.  DS. DIST  EELEC. ETON. ICELL. 1 STOF . 


+  IDS,  THEAT,  ALF'HS) 

IF ( I STOP. EQ. 1 )  ICELS=ICELL 
IF  (ISTOF'.EQ.  1 )  GO  TO  600 
XO=XN 

THETO=THETN 

500  CONT I NUE 
600  CONTINUE 

HE ( I CELL ) =HE ( ICELL ) +ALF HS 

AMASS (ICELL) =AMASS ( I  CELL) +ALPHS*6 . 635E-24 
3000  CONTINUE 

GO  TO  2000 
250  CONTINUE 

ALPHS=REAC (I  TONE) +ALPHA ( I  ZONE) 

HE  ( I  ZONE)  =HE  <  I  ZONE)  +AL.PHS 

AMASS  ( I  ZONE )  =  AMASS  ^  I  ZONE )  •+•  ALF'HS ♦  6 . 635E— 24 
DENSE  =  AMASS ( I  ZONE / /V ( I  ZONE) 

FI  =  1.  /  (l.  +  (3.2.-.B63*AL0e<DENSE>  )  /TE(IZONE)  ) 
E I  ON  ( I Z  ONE )  =3500 *  ALPHS * F I 
EELEC ( I  ZONE ) =3500 . RALPHS* ( 3 . -F 1 ) 

THEAT ( I  ZONE) =3.5* 1 . 602E-6*ALF'HS 
2000  CONTINUE 


3600=C 


3610=C******  COMPUTE  NEW  TE,  TI 


3620=C 

3630= 

3640= 

3650= 

3660= 

3670= 

3680= 

3690= 


DO  700  I ZONE= 1 . NSPACE 

ET0TE=3. /2. *TE ( I  ZONE ) *ELEC(IZONE> 

ETOTE---ETOTE+EELEC  ( I  ZONE) 

TE ( I  ZONE) =ET0TE*2. /3. /ELEC (I  ZONE) 
ET0TI.=3.  /2.  *TI  (I  ZONE)  *H(IZONE) 
ETOTI=ETOTI+EION( I  ZONE) 

TI < I  ZONE) =ETOTI *2. /3. /H(IZONE) 


3700= 
37 1 0= 
3720= 
3730= 
3740= 
3750= 
3760= 


700  CONTINUE 
RETURN 
END 

SUBROUTINE  \YT3 < XO, XN.  DIST. THETO, THETN, INOW.R) 
DIMENSION  R(101) 

RR=R  ( INOW+1  > 

XN=RR 


3770=  ANG=ACOS (XO/RR) 


IAL 
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37 oO~  V7  E  MP=  F  P  <  I N  ( r.N5 ) 

77-7  !'=•  7 HETN-THETQ— AN3 

38  »>=  DU'l  =V TE-'iF' 

3810=  IN0W=IN0W+1 

3820=  RETURN 

3830=  END 

394 SUBROUTINE  XYT4  ( XO,  XN,  DIST  THETO.  THETN.  INOW,  R) 

3950=  D I MENS I ON  R (101) 

3960=  CALL  BETA (THETO. A) 

3970-  CALL  SETS (XO. THETO.  EO 

3830-  F'R=R i INOW) 

3890=  RP 1 =R ( I NQW+ 1 ) 

3900=  CALL  GETC  (FF. .  XO  THE1U.O 

3910=  CALL  GETC'RRi , XO. THETO. Cl ) 

3920=  PISC=B**2-4. *A*C 

3930=  DISCl=B**.2-4.  *A*C1 

3940=  IF (DISC. BE. 0. )  CALL  XYT4A (THETO,  RR.  XO.  B,  DISC, A, XU. THETN. DI3 


3950= 

3960= 

,  N, 
3970= 
3980= 
3990= 
4000= 
CELL) 
4010= 
4020= 
4030= 

4  040= 

4050= 

4060= 

4070= 

4080=0 

4  0~'0= 

4  J  00= 

4  1 1  0= 
4120= 

I FLAG, 
4130= 
4140= 
4150= 

4  160= 
4170= 

4 1 80= 

4 1 90= 
4200= 
4210= 
4220= 
4230= 
4240= 
4250= 
4260= 
427r,= 
4280= 


(-  O,  I  DUN) 

IF  (DISC.  LT.  O.  )  CALL  X  YT4B  <  INOW ,  RR  1 .  THETO .  XO,  A .  D I  SC  i ,  B .  i'HETN 

*■  DIET,  0,  I  DUN ) 

RETURN 

END 

SUBROUTINE  XYT4A (THETO, RR.  XO,  B.  DISC.  A.  XN.  THETN, DIST, IFLAS. I 

F 1=3. 141592654 

YTEMF--  (-B-SORT  (DISC)  >  /2.  /A 

XTEMF‘=  ( YTEMP— XO+TAN  (PI— THETO)  >  /TAN  (THETO) 

DIST=SQF:T  <  (XO-XTEMP) **2*YTEMP**2> 

XN=RR: 

THETN=THETO- AS I N ( YTEMF / RR ) 

IC£LL=IN0W-1 

Il-.!0W=JN0W 

I F  ( I FL  AG .  EO.  1  )  I  MOW*-- 1  NON- 1 

RETURN 

END 

SUBROUTINE  XYT4B (IN0W.RR1 ,  THETO,  XO,  A.  DISCI . B, THETN. XN. DIST, 

i-  I  CELL) 

PI =3. 141592654 

YTEMP=(~B+SQRT  (DISCI) ) / 2 .  / A 

XTEMF- (YTEMP-XO* TAN (FI-THE TO) ) /TAN (THETO) 

DIST “SORT ( (XO— XTEMF ) **2+YTEMF**2> 

XM=RR1 

THETN=THETO-ASIN  (YTEMP/RR1 ) 

ICELL=INOW 

IN0W=IN0W+1 

IF ( IFLAG. EQ. 1 )  INOW=IMOW-1 

RETURN 

END 

SUBROUTINE  XYT5(X0, XN. DIST, THETN. IWOW.R) 

DIMENSION  R(101) 

PI =3. 14159265^ 

XN=R ( INOW) 


112 


r 


4300= 
431  0= 
4 330- 
4330= 
4340= 
4350= 
4360= 
437 0= 
4380= 
4390= 
4  4  00= 
4410= 
4  420= 
4430= 
4440= 
4450= 
4460= 
4470= 
4  480= 
4490= 
4500= 
4  5  1 0= 
4520= 
4530= 
4540= 
4550= 
4560= 
4570= 
4580= 
4590= 
4600= 
4610= 
4620= 
4630= 
4640= 
4650= 
4660= 
4670= 
4680= 
-3  690= 
4700= 
4710= 
4720= 
4730= 
4740= 
4750= 
4760= 
4770= 
4730= 
4790= 
4800= 
4810- 
4820= 
4  870= 


L 


*•'  .-FTW=  F  I 
I N0W= I NOW 
D I SI =ABS  <  XN— XO) 

J  F ( R ( I NOW > .  LE. 0 .  )  THE  TN=0 . 

RETURN 

END 

SUBROUT I NE  BETA (  THETO . A ) 

A=1 . +1 . / (TAN (THETO) >  **2 

RETURN 

END 

SUBROUTINE  GETS ( XO. THETO. G > 

PI =3. 141592654 

ANUM=-2 .  *  XQ  it  TAN  ( P I  -THETO ) 

DEN= (TAN (THETO) ) **2 

B=ANUM/DEN 

RETURN 

END 

SUBROUTINE  GETC  <R. XO. THETO, C) 

PI =3. 141592654 

ANUM=  ( XQJfcTAN <F' I  —THETO )  )  **2 

DEN= (TAN (THETO) ) **2 

C=-l . *R**2+ANUM/DEN 

RETURN 

END 

SUBROUTINE  XYT6 ( I HOW. R. XN.  THETN. DI3T. I CELL, XO) 
DIMENSION  R(101) 

XN=R ( INOW+1 ) 

THETN=0. 

I CELL* I NOW 
I NOW= INOW+1 
DIST=XN-XO 
RETURN 
END 

SUBROUTINE  XYT7(X0. I NOW,  THETO.  R, DIET. ICELL. XN, THETN) 
DIMENSION  R(101) 

F'I=3.  14  1592654 
CALL  BETA (THETO. A) 

CALL  GETB(XG. THETO. B) 

RR=R( INOW+1) 

CALL  GETC (RR.XO. THETO. C) 

YTEMF- (-B+3QPT (B* *2-4 . *A*C> ) /2. /A 
XTEMP=  <  YTEMF'-XO*TAN  (PI-THETO)  )  /TAN  <  THETO) 
THETN=THETO-ASIN  (YTEMF'/RR) 

XN=RR 

ICELL=INOW 
I NOW= INOW+1 

DIST=SQRT < (XO-XTEMF) **2+YTEMP**2) 

RETURN 

END 

SUBROUTINE  XYT8  *  I NOW, THETO. XO.  THETN. XN.DIST.R) 
DIMENSION  R(101) 

CALL  GET  A (THETO, A) 

CALL  GETE <X0. THETO. B> 

F:R=R  ( I  NOW— 1  ) 

RP 1 =R ( I NOW ) 
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r 


■IR'l'.  - 
4850= 
4860= 
4370= 
4  330= 
T,  1, 
4890= 
4900= 
.  XN. 

4  910= 
4920~ 
4930= 
4940= 
4950= 
4960= 
4970= 
4980= 
4990= 
5000= 
5010= 
5020= 
5070= 
5040= 
5050= 
5 060= 
5070= 
5030= 
5090= 
5100= 
5 1 1 0= 
5120= 
5130= 
5140= 
5150= 
5160= 
5170= 
5180= 
5190= 
5200= 
OP, 
5210= 
5220= 
5230= 
5240= 
5250= 
5260= 
5270= 
5280= 
5290= 
5300= 
5310= 
5320 
533<'>= 
5340= 
5350= 


’  3r"  fnF'.  ,-l»  THF  T 0.  C) 

CALL  C5ETC  (RF  1  -  XO.  ThETO.Ci  > 

D1 SC=B**2-4. *  A*C 
DI5Cl=B**2-4.  *A*C1 

IF  CDISC.GH.O. >  CALL  vYT4A (THSTO,  RFL  XO.  3. DISC, A. >N. THETN. DIS 
+  I CELL) 

IF  (DISC.LT.O.  )  CmLL  XVT4B  ( I  NOW.  PR  1 ,  THETO.  XO.  A.  DISCI .  B.  THETN 

+  '0157 .1.1  CELL ) 

RETURN 

END 

SUBROUTINE  XYT9 ( TNOW. XO. THETO.  R,  1 CELL.  D I ST. XN. THETN' 
DIMENSION  ROOD 
THETN=THETO 
XN=R ( INOW-1 ) 

DIST=ABS (XN-XO) 

I CELL= INOW-1 

I WOW= INOW-1 

IF  (  XN.  EG!.  0.  )THETN=0. 

RETURN 

END 

SUBPOUT I ME  F I HDDS i TEMPE .  DENSE . DS . D I ST , IDS ) 

RATIO®. 213/PENSE 

ANUME=TEMPE** 1 . 5*0. 036+RATI0 

DEME=1 . +0. 17*AL05(TE!-1PE*SQRT  (RATIO) ) 

DE=ANUME/DENE 

ANUMI=10.65*RATI0 

DENI  =  1 .  +0.75* ALOE  (SORT  (  TEi-IFE* RATIO)  ) 

DI=ANIJHI/DENI 
IF  (TEMPE.  L-E.  20.  )  D=DE 
IF (TEMPE. BT. 20. )  D=DI 
D=D/50. 

IDS=INT (DI5T/D) 

IF(IDS.LE.l)  IDS=1 

DS=DIST/IDS 

RETURN 

END 

SUBROUT  I TJE  EL05S  ( DENSE ,  TEMPE  .  U .  DS .  D I  ST ,  EELEC ,  E I  ON .  I  CELL  .1ST 

+  I ALF  H  .IDS.  THE  AT ,  ALF’HS ) 

D I MENS I ON  THE AT (101) 

D I MENS I ON  EELEC ( 1 O 1 ) . E I ON (1 0 1 ) 

SF'E  (D,  U„  T)  =—23. 2 /D* SORT  (U)  /T**1.5*  U.--.  17*AL0G  (T*SQPT  (D)  )  > 
SPI (D, U, T) =—.047 /D/Ut <  3  .  +0. 075* ALOG (SORT (T*D*U>  )  ) 

RAT  1 0=. 213/DENSE 

IF (RATIO. EQ.O. . OR. TEMPE. EG. O. )  STOP 
DO  100  1=1, IDS 
FRACE=SPE (RATIO, U, TEMPE) 

FF.'ACI  =5F'I  (RATIO,  U,  TEMPE) 

U=U+DS* (FRACE+FRACI) 

UTHEF:M=3.  /2.  ''3500.  t TEMPE 
IF ( U . LE . UT  HEPM '  IST0F  =  1 
IF (U. LE. UTHEFM)  GO  TO  200 
EE=3500. *D5*AB5 (FRACE) *ALPHS 
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...i-c'O” 

5390= 
ALPHS 
--  4,30= 

5-;  •  0  ~ 
5420= 


440=  3 


..  54  <K  ~ 


MAE  = , 

St.  90= 
5700= 


t  "50'  *  *AM  (.-RAM )  *  Lr'H2 

F  9  93  ( ]  CELL  >  =[>!  SC  ( ICE!  L)  "-EL 
t  I  Oh  <  I CF  L  •  =E  I  'JW  (  I  CELL  5  *  3  r 

THEr-il  ( I  CELL )  =THEAT  t  I  CELL  >  •+•_■»  5*DS*  ABb  ( r  RACE-'T-RAC  ]  i  f-j. ,  c'.CC-o* 

I F  < 1 ALF  H . EQ . 1 . AND . E 1 . GE . EE  >  GO  TO  300 

IOC  CONTINUE 
200  CONTINUE 
RETURN 
300  15409=1 

E I  ON  1 1  CELL )  =EJ  ON  <  I CEI  L  >  -<-U  *  3  5  >0  -  *  ALr’HS 

7  HEAT  (  T CELL )  =1  HEAT  ( I  CELL)  t  3.  5*  1  -  £  02E-6:f  U RALPHS 

RETURN 

END 

5  -=0  '~TUE  HTFLX (R.7E.TJ.7R  3LEC  V. H. HE. PUSH. A MASS. NSPAUt, H 
*  MEMASG. PM ASS. DT, PR) 

DIMENSION  F:  (101 5  „  TEC  101  .TR(101)  ,TI  <301)  .ELEC  (101)  , 

V < 1 o 1 5  .  PR (101).  AMASS (101) . H ( 1 0 1 ) , HE (101), PUSH  ( 1 C 1 ) 

F  E.  v.  LOGLAK .  HF  .HI .  K .  KP .  MU .  NUB .  NuC .  URGE .  LAMROS .  NT  2 OVA ,  NEP .  MG 


HEAT  CONDUCTION  13  CONSIDERED  USING  COEFFICIENTS  TAKEN  FROM 
FRALEY.  ET  AL.  FHYS  FL.  VOL  17.  NO  2,  FEB  7-..  USING  FLUX-CONSE 


RVATIVE 

=- .  —  -  — n  r  t: 


5^.0  «C  ELECTRON  MASS 
5770=  mi  . -iiL-r: 

57S0«C  BOLTZMANNS  CONSTANT  IN  ERGS/ KEV 
5790=  K=1 . 39E-13* 1  lalE+7 

5300=0  ELECTRON  CHARGE 
57  1  0~  E=4 .  BE— 1  0 


5020=0  FINE  STRUCTURE  CONSTANT 
5330=  ALF-1/137.04 

3940=  F 1=3. 14159 

I3S5'>=C  SPEED  OF  LIGHT 
5830=  C~2. 9979E10 

5 370=0  STEF'HAN-BOLTZMANN  SIGMA  IN  ERG3/CM2/SEC/KEV4 
5330=  SIGMA  =  5. 67E-5* ( 1 . 161E+7) **4 

3390=  FLUXM=0. 0 

5900=  RFLUXM  =  0.0 

5910=0  COLLECT  ALL  OF  THE  CONSTANTS  IN  FRONT  OF  THE  DG-LOOP 
5920=0  A  AND  B  ARE  USED  IN  LOG  LAMBDA  =  LOGLAM  (EON  R4C) 

5930=C  D  IS  USED  IN  ELECTRONIC  HEAT  DIFFUSION  COEFFICIENT  .=  KP  (EON 
ETC) 

5940=C  F  13  USED  IN  ELECTRON- I ON  RELAXATION  RATE  =  NU  (EDN  P4A) 
5950=C  HI  AND  H2  ARE  USED  IN  ELECTRON-RADIATION  RELAXATION  RATE 
5960-C  =AERB  +  AERC 

5970 =C  FLUXP  AND  FLUXM  ARE  EACH  AN  ELECTRON  HEAT  CONDUCTION  FLUX  TIM 

ES  AREA 

59E0--C  RFL.UXP  AND  RFLUXM  ARE  SIMILAR  FLUXES  TIMES  AREA  FOR  RADIATION 
5990=  A=(l . 5/E**3> *»  l . 5/SQRT(PI ) 

9000=  B=. 5/ALF7C*SQRT  (3.  /ME) 

6010=  D=20.  *  (  (2.  /PI)  **1.5)  *!:V  (SOFT  (ME)*  (E**4)  ) 
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4020=  F  =  5.  0tSC-RT(2.  +  F  T  *ME)  *E*  *-1/2. 

4030-  HDARC=6.  soE-27'2.  'FT  *C 

60m  0=  MC2=HE*C.tC 

6050=  RO=E  •*  E/MC2 

6060=  H 1  =32.  '7.  *SQRT  <  <2. /PI /ME)  >  *  <E»-E*E*E)  *1-  *R0/KEAF:C 

6070=  H2=  (128.  /3.  )  *PI  *  <E  'MCI.)  *  (E/MC2)  +5IGMA*R0 

6080=  DO  ICO  N=1 , NSFACE 

60?0=C  ELECTRON- I ON  EXCHANGE  AND  ELECTRON  HEAT  CONDUCTION 
6100=  RF  -F;  (N) 

6 1 1 0=  2 2=ELEC  ( N ) / (H  < N >  +HE  < N )  ■♦•PUSH  ( N )  ) 

6120=  IF  (N.LT.NSPACE)  GO  TO  10 

6130=  TP=TE  CN) 

6 1 40=  NEF=EL£C ( N ) / V (N) 

6150=  ZP=22 

6160=  GO  TO  20 

6170=  10  TF'=.  5*  (TE  (N+T  >  +TE  (N)  ) 

6180=  NEP=.5*  (ELEC(NJ-1)  /V (N+l )  +ELEC (N)  /V (N)  ) 

4190=  Z  1  =El.EC  ( N+ 1)  /  ( H  (N+ 1 )  +HE  ( M+ 1 )  +PUSH  < M+ 1 )  ) 

6200 ~  ZP=.  5.f  (Z1+Z2) 

6210=  20  EOT  =  .  43*ZP/ <3.  44+ZF'+.26>:A!-0E  <7.p'‘ ) 

6220=  SQTE=SG)RT(TP) 

-.23  '  =  LOSLAM=Ai  OG  ( BOTE *  *  3/NEP /  <ZP+-R*SQTE)  > 

6240=  IF  (L03LAM.LT. i  LQGi_AM=l. 

6250=  KP«D*  (TP**2.5)  *EDT/ (  ZP.fLOGLAM) 

6300  MI  =  'HMASStH  (N>  +HEMASS*HE  (N)  -PMASS*PU«H  <N>  )  / 

63 1  - )— ■  /.  ( H  ( M )  =HE  <  H )  +PU3K  ( N )  ) 

6320=  NU=F*ELEC  ( i  i )  *FLEC  IN)  *LOGLAM  f  (AMASS  (N)  *V  (N)  %TE  <N)  ol.  5 

6330=  IF (N. 5F.NSP ACE-2)  GO  T0  22 


6340= 

DTDRP 

6350= 

GO  TO 

6360= 

22  DTDF'F= 

6370= 

23  FLUXP- 

6330= 

i'E  (N) 

6390= 

1  -t-Nl 

l  400= 

TT CM) = 

64 10= 

FLUXM 

6420=C 

RAD I AT I C 

6430= 

IF  (N. 

6440)= 

TRP  = 

6450= 

TRP3  = 

6460= 

DTDRP 

6470= 

rEF  = 

4  4  30- 

~"CT  pi*  T  r? 

64-'0=  '-'HOF  2  —  (AMASS  (N)  /V  (N)  )  *#2 

6300-  9p-R  •;;•!)  *R(M) 

4710~  GO  TO  33 

6520=  32  TRr  =  (TR (N+l) =TT (N) ) /2. 

6530=  TFF3  =  TEP‘TRP*TRF 

65^0=  2TDPF  =  <TR<N=2>  -TR(N)  )/ (R  (N+2) (N) 

6550=  TEC  =  (TE  (M-t-1 ) +TE  (N)  ) /2. 

6560=  PRINT*, TEP 

6570=  TEP72  =  TEF+*3.5 


i>°  c0° 

ay 


0  i  N*- 1 )  +AMA5S  -U)  /V  (N)  )  *0  a. 
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oaM  \j  t  i  —  .  *  i  R  1 1 )  /  T  E  t  -i  i 

6650=  G=F  I  tPI/4  .  -  (PI*F  1/4.  -1  >  *EXP  l-Al) 


6670= 

TR4  =  ( T R  ( N >  *  2 ) 

*+2 

6680= 

TRT  =  TR4/TF;(N> 

6700  = 

A”rF'H , /SCR 7  TE ( 

N>  )  *  (ELEC  i 

•J)  »HL£ 

C  ( ! ! ;  V  !■’;  >  /V 

•£:  T  i  0  ~ 

AE?;.C=H2t NER *  TF4 

c ‘  u i  -  i.5vf.,M: 

E,  "7  “»■  (m/ 

C'  1E  =  CNU I  *  72 

.t”40= 

■  SF=,-iFRE  +AEPC 

6751  ■= 

C  N  U  F  i  -t  3  I G  M  A  /  C  t  V  ( !  4 )  /  A  M  r  i  S  i 

CN)  .1  TF, 

-r 

te  (H)  =  to  an  - 

AER  y (TE (N) 

-  tp  (  m 

!  iDT/CNL'E 

677':= 

TF  ( N  )  =  TF;  ( N  >  + 

Ar  R*  (TE  ( N > 

— TF  (N) 

)  IDT  'Ct  UR 

6780=  '  + (RFLUXP-PFLUXM) *4 . *PI *DT  ( C (N) *CNUR) 

6790=  RFLUXM  =  RFLUXP 

6510=  100  CONTINUE 

6320-  RETURN 

^C3 1-  I', 

6680=  SUBROUT I NE OUTPUT  (R.  CYCLE,  TIME.  D7.  DREl  5.  NRJEL.  N5PACE.  VEL  ,  RHO 

Z . PM AS 3 

6690=  /.,  ZP’UEiH,  ELEC.  TE,  TI ,  TP.  RHAC.  TREAC ,  THEAT  .  FLU .  RELPHD.  H.  HE,  PUSH, 

DRVENG > 

6700=  D I  MENS  7  ON  VEL  (101),  PELF; HD  ( 1 0 1 )  .  H  ( 1  0 1  .  HE  ( 1 0 1 )  ,  PUSH  ( 3  0 1  )  .  ELE 

C(101)  . 

-.710-  TE  (101  >  .  TI  (lOl)  .  TP  (101),  READ  (101) .  TPEAC  <  1  01 )  .  THEAT  ( 101 )  .  PL 

U(  101 )  . 

6~C0-.  DP  VEND  ( 101  )  .  R  <  3  0 1 ) 

6730=  INTEGER  CYCLE 

6740=  TPL0T=TIME*1.E9 

6750=  DTF'L0T=DT*1 .  E12 

6760=  '  E I PL0T=DREN6* 1 .  E- 1 0 

6770=  E0PL.0T=0, 

6780=  DO  2  M-=J  .  MFL'EL 

6790=  2  EOPLOT=EGPLOT+REAC(N> *2. 816E-15 

6800=C  EOF'l  QT  IS  LJ  RELEASED  Bv  DT  BURN 

6810=  PRINT  101.  CYCLE .  TF'LQT ,  DTPLOT ,  E I  PLOT ,  EOPLOT 

6820=  PRINT  102 

6830=  PRINT  104 

6840=  DO  10  N=1  .  NFUEL 

6050=  RPLOT=F;  ( N )  *  1 0 . 

6860=  VpUn=VEL(N>  * l.E-8 

6870=  OUTFL=TREAC (N) *2. S16E-12 

6880=  HTPLT =TMEAT (N) *  3 . E-7 

6890=  10  PRINT  103, N,  RPLOT,  VFLOT,  REAC  (N>  .  REl.RHO  (N)  .  TE  <N>  .  TI  <N>  ,  TF;  <N> 

6900=  V.OUTF'L  .  HTPLT,  FLU  (N) 

6910=  IF  (NFUEL,  LE.NSF’ACE)  GO  TO  200 

6920=  PRINT  105 

6930=  NFF-NFUEL+1 

6940=  DO  20  N=NFF‘ ,  NSPACE 

6950=  F'F'LOT =R  ( N )  ♦  1 0 . 

6960  -  VF'LOT =VEL  (N)  *  1  .  E-B 
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DU  TPLOT=r*PVENP  i  N  *  1  .  E-? 

-  HTPLT -THi  "iT  H)  i  '  .  F 

»9?0=  20  PRINT  i  03,  N,  RPLOT .  VPLOT ,  REAC  •.  (i )  .  RELRHO  ( N )  ,  IE  «N>  .  T  J 


•  ‘Oc= 
"01 0« 
7  20= 
~030- 

?050= 
1 1>0= 
-r<  >70- 
’  /  / 


7.0UTPL0T.  HTF'LT,  FLU  (N> 

1 0 1  FORMAT  < ’  1  CYCLE  = ’ .  I A . 3  X .  ’ T I ME  = ' . FB . 4 . ’ ME  DT 

7.F8.4.  ’FS  INPUT  =  '.F3.4.’f  7  OUTPUT  =  7  .  ^3.  3.  7  f  J’  . /, /) 

1  02  FORMAT  1 7  >: .  ’  Rh t  >  l  US  0  >:  7  VEL  0'.'  I T  \  .  T  .  7 !  EACT I  DUS  7  . 

1  T  a  .  7  DEN  -  I T  Y  '  .:.X.  7  Ei ..  1EHP'  .4/..  •’  I  DM  TE  MF'  ’  .  <1  X .  7  F  L>  TFMC  ’  . 

2  5  X .  ’OLD  ,  UT7 .*X.  ’  HEAT  IMS 7  .  EY.  7  NEU'i  RON’  ./.3>  7  1 7 . -X . 7  (  MM  > 

3  EX.  ’  (MM/MS*  ’  .4X.  'THIS  I'T'  .  3X  ’  (X  N~Rt  Af  )  ’  .  -•>  '  (KE>'>  7  . 
47X.  ’  <KEV)  7  ,7X.  7  (LEV'  •'  ,5X,  7  (JOULES)  7  -2V.  ’  (JOULES)  7  .  5X.  'FI 


7080= 

5) 

7090= 

103  FORMAT (IX. 

I5.2F12.6.  1  PE  12.4.  OF1 . 6F 1 2 . 3 ,  1PE12 

7  1 00= 

1 04  FORMAT  ( 7 

FUEL 7 , //> 

^1  10- 

105  FORMAT  (// 

.  ’  PUSHER  ’  .  90X .  7  I  IMF'UT '  .  /  /  ) 

"  1  20= 

200  CONTINUE 

"  1  30= 

RETURN 

7 1  40= 

END 

,  EN 
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APPENDIX  B 
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NFUEL 

NSPACE 

PMASS 

PUSH(I) 

PUSHD 

R  ( I ) 

TE  (I) 

THEAT(I) 

TI  (I) 

TIME 

TMASS 

TR  (I ) 

TREAC(I) 

VCD 

VEL(I) 

ZPUSH 


Variable 

DELV(I) 

DR 

NSSP1 


total  number  of  fuel  zones 
Total  number  of  zones 
Pusher  mass  in  grams 

Total  number  of  pusher  atoms  in  zone  i 
Pusher  material  density  in  grams 
per  cubic  centimeter 

Radial  distance  of  zone  i  in  centimeters 
Electron  species  temperature  in  kev  in 
zone  i 

Total  heating  in  ergs  of  zone  i 

Ion  species  temperature  in  kev  in  zone  i 

Time  in  seconds 

Tritium  mass  in  grams 

Radiation  temperature  in  zone  i 

Total  number  of  reactions  in  zone  i 

Volume  in  cubic  centimeters  of  zone  i 

Velocity  of  cell  wall  i  in  centimeters 

per  second 

Average  atomic  number  of  pusher  material 
Table  10 

Subroutine  GDATA  Variables 

Definition  or  Description 

Change  in  volume  in  cubic  centimeters 
Zone  width  in  centimeters 
NSPACE  +  1 
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PIE 
PR  (I) 

RM3 

RP3 


Variable 

A 

ACC  (I) 

AD 

AN 

C 

DIV 

DTI 

DTMAX 

EFFTE 

FRACT 

PE  (I) 

PI(I) 


IT 

Radiation  pressure  in  ergs  per  cubic 
centimeter 

Inside  zone  radius  in  centimeters  cubed 
Outside  zone  radius  in  centimeters  cubed 

Table  11 

Subroutine  HYDRO  Variables 
Definition  or  Description 

Constant  =  407'c  for  radiation  pressure 
calculation 

Acceleration  of  zone  wall  i  in  centimeters 
per  second 

Denominator  for  acceleration  computation 

Numerator  for  acceleration  computation 

Speed  of  sound 

Divergence  of  velocity 

Interval  time  step 

Maximum  time  step 

Effective  electron  temperature 

Fractional  value  for  time  step  calculation 

Pressure  in  ergs  per  cubic  centimeter  due 

to  electron  species 

Pressure  in  ergs  per  cubic  centimeter  due  to 
ion  species 
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PIE 

PR(I) 

PT(I) 

Q  ( I) 

RELRHO (I ) 

RHO(I) 

RHOOLD 

SIGMA 

VEL(I) 

VOLD 


Variable 
I  BURN 

SIGVDT 
TM13 
TM23 
T START 


IT 

Pressure  in  ergs  per  cubic  centimeter  due 
to  radiation 

Total  pressure  in  ergs  per  cubic  centimeter 
Artificial  viscosity  in  ergs  per  cubic 
cent imeter 

Relative  density,  current  density  over 
ambient  density 

Density  in  grams  per  cubic  centimeter 
Density  during  previous  iteration 
Stephan-Boltzmann  constant  in  ergs  per 
square  centimeter  per  second  per  kev^ 

Cell  wall  velocity  in  centimeters  per 
second 

Volume  during  previous  iteration 
Table  12 

Subroutine  TBURN  Variables 

Definition  or  Description 

Logical  flag  used  to  carry  burn  status 
to  other  subroutines 

<av>DT 

Cubic  root  of  ion  temperature 
Ion  temperature  to  the  2/3  power 
GO/NO  GO  temperature  for  thermonuclear 
burn 
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Table  13 

Subroutine  NUHEAT  Variables 


I 


Variable 

ACTION 
ALPHA (I) 


ANO(I) 

AVGEN(I) 

BNO(I) 

EE 

EX  (I) 

S 

SIGMA 

TOTINT 


Definition  or  Description 

Total  number  of  interactions 
Equivalent  number  of  3.5  Mev  alpha  particles 
computed  from  energy  transfer  from  neutron- 
ion  collisions 

Ion  number  density  in  zone  i 
Average  energy  lost  by  a  neutron  during 
collision  with  an  ion 
Total  ion  number  in  zone  i 
Probability  neutron  penetrates  zone  i 
Probability  neutron  penetrates  zone  i+1 
Total  number  of  neutrons 
Neutron-ion  cross  section 
Total  number  of  neutrons 


Subprogram  ALPHA1  is  discussed  by  subroutines. 


Table  14 

Subroutine  ALPHA1  Variables 
1.  Subprogram  ALPHA1 

Variable  Definition  or  Description 


ALPHS 

DENSE 

EELEC(I) 


Number  of  new  alpha  particles 

Zone  density  in  grams  per  cubic  centimeter 

Energy  transferred  to  electrons  in  zone  i 


9 
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E I  ON  Cl) 
ETOTI 
ETOTE 
FI 

IALPH 

ICELL 

ICELS 
I  LOOP 
INOW 

ITHET 
I  STOP 

I  ZONE 

TEMPE 

THET 

THETA (I) 

THETN 

THETO 

U 

XN 

XO 


Energy  transferred  to  ions  in  zone  i 

Total  energy  of  ions 

Total  energy  of  electrons 

Fraction  of  energy  deposited  in  ions 

Logical  flag  for  energy  deposition  position 

Zone  number  into  which  alpha  particle  is  about 

to  decelerate 

Vestigial  statement 

Simple  loop  parameter 

Current  alpha  particle  position  on  radial 
mesh 

Index  for  angular  mesh 

Logical  flag  to  signal  alpha  particle 

thermalized 

Radial  zone  number 

Electron  temperature 

Particular  angle  of  angular  grid 

Angles  of  angular  grid 

New  theta  angle 

Old  theta  angle 

Remaining  fraction  of  alpha  energy 
New  x  coordinate 
Old  x  coordinate 


124 


2.  Subroutine  XTY3 


Variable 

ANG 

DIST 

RR 

YTEMP 

3.  Subroutine 
Variable 

DISC 

DISCI 

RR1 

4.  Subroutine 
Variable 
XTEMP 


Definition  or  Description 

Angle  measured  from  microsphere  center 
from  birth  point  to  exit  position 
Distance  traveled  within  the  zone 
Radius  of  zone  being  entered 
Distance  traveled  within  the  zone  in  the 
y  direction 

XYT4 

Definition  or  Description 

Discriminant  of  the  solution  of  the 
equation  of  the  line  of  travel  with  the 
circle  representing  the  next  inner  zone 
Discriminant  of  the  solution  of  the  equation 
of  the  line  of  travel  with  the  circle 
representing  the  next  outer  zone 
Radius  to  the  outside  of  the  zone  being 
entered 

XYT4A 

Definition  or  Description 

Distance  traveled  within  the  zone  in  the 
x  direction 


125 


5.  Subroutine  XYT4B 


No  New  Variables 


6.  Subroutine  XYT5 
New  New  Variables 

7.  Subroutine  GETA 

Variable  Definition  or  Description 

A  Coefficient  of  the  solution  of  the 

intersection  of  a  line  with  a  circle 

8.  Subroutine  GETB 

Definition  or  Description 

Numerator  for  coefficient  calculation 
Coefficient  of  the  solution  of  the 
intersection  of  a  line  with  a  circle 
Denominator  for  coefficient  calculation 


Variable 

ANUM 

B 

DFN 


9.  Subroutine  GETC 

Variable  Definition  or  Description 

C  Coefficients  of  the  solution  of  the 

intersection  of  a  line  with  a  circle 
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10.  Subroutine  XYT6 


No  New  Variables 

11.  Subroutine  XYT7 
No  New  Variables 

12.  Subroutine  XYT8 
No  New  Variables 


13.  Subroutine  XYT9 
No  New  Variables 


14.  Subroutine  FINDDS 


Variable 


Definition  or  Description 


ANUME 

ANUMI 

D 

DE 


DENE 


Numerator  for  range  calculation  if 
electrons  only  decelerate  alpha  particle 
Numerator  for  range  calculation  if  ions 
only  decelerate  alpha  particle 
Intermediate  variable  supporting  integration 
increment  calculation 

Range  of  alpha  particle  in  zone  if  zone 

goes  to  infinity  and  only  electrons  decelerate 

it 

Denominator  of  range  calculation  if  electrons 
only  decelerate  alpha  particle 
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DENI 

DI 


DS 

IDS 

RATIO 

15.  Subroutine 

Variable 

EE 

El 

FRACE 

FRACI 

SPE 

SPI 

UTHERM 


Denominator  of  range  calculation  if  ions 
only  decelerate  alpha  particle 
Range  of  alpha  particle  in  zone  if  zone 
goes  to  infinity  and  only  ions  decelerate 
it 

Integration  increment 

Number  of  integration  increments  within 
current  zone 

Relative  ratio  of  deuter ium-tritium 
compared  to  ambient 


ELOSS 

Definition  or  Description 

Energy  to  electrons  in  kev 
Energy  to  ions  in  kev 
Electron  term  for  energy  deposition 
Ion  term  for  energy  deposition 
Electron  term  for  energy  deposition 
Ion  term  for  energy  deposition 
Energy  of  thermal  ions 
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Table  15 


Variable 

A 

AER 

AERB 

AERC 

ALP 

A1 

B 

CNUE 

CNUI 

CNUR 

D 

DTDRP 

E 

EDT 

F 

FLUXM 


Subroutine  HTFLX  Variables 
Definition  or  Description 

Constant  for  Coulomb  logarithm  computation 
Radiation-electron  energy  transfer  coupling 
coefficient 

Bremsstrahlung  coupling  coefficient 
Compton  scatter  coupling  coefficient 
Fine  structure  constant 
Constant  for  bremsstrahlung  Gaunt 
factor  calculation 

Constant  for  Coulomb  logarithm  computation 

Specific  heat  for  electron  species 

Specific  heat  for  ion  species 

Radiation  specific  heat 

Constant  for  coefficient  of  thermal 

conductivity  calculation 

Spatial  derivative  of  temperature 

Electronic  charge  in  esu 

Term  e6^,  for  coefficient  of  thermal 

conductivity  calculation 

Constant  for  equilibration  coefficient 

computation 

Heat  flux  at  zone  inner  radius 
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FLUXP 


Heat  flux  at  zone  outer  radius 


G 

HI 

H2 

KP 

KROS 

LAMPOS 

LOG LAM 

MI 

NEP 

NU 

RFLUXM 
R FLUXP 
RHOP2 
RP 

SQTE 

TEP 

TEP72 

TP 

TRP 

TRP3 

TR3 


Bremsstrahlung  Gaunt  factor 

Constant  for  bremsstrahlung  coupling 

coefficient  calculation 

Constant  for  Compton  scatter  coupling 

coefficient  calculation 

Coefficient  of  thermal  conductivity 

Diffusion  coefficient  for  radiation 

Mean  free  path  of  photons 

Coulomb  logarithm 

Average  ion  mass 

Electron  density 

Equilibration  frequency  for  electron-ion 
energy  transfer 

Radiation  energy  flux  at  inner  zone  radius 
Radiation  energy  flux  at  outer  zone  radius 
Square  of  zone  density 

Cell  wall  radius  from  center  in  centimeters 
of  zone  in  loop 

Square  root  of  electron  temperature 
Electron  temperature  in  kev  of  zone  in  loop 
Electron  temperature  to  7/2  power 
Average  electron  temperature  over  two  zones 
Radiation  temperature  in  kev  of  zone  in  loop 
Cube  of  radiation  temperature 
Cube  of  radiation  temperature 
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TR4 

Radiation  temperature 

to 

the  fourth 

power 

ZP 

Net  electronic 

charge 

averaged 

over 

zones 

i  and  i+1 

Z1 

Net  electronic 

charge 

in 

2  one 

i+1 

12 

Net  electronic 

charge 

in 

zone 

i 

Table  16 

Subroutine  OUTPUT  Variables 


Variable 

Definition  or  Description 

DTPLOT 

Time  step  in  picoseconds 

E1PL0T 

Energy  in  kilojoules  input  by 

dr  iver 

EOPLOT 

Energy  in  kilojoules  released 

by  thermo- 

nuclear  burn 

HTPLT 

Bootstrap  heating  in  kilojoules 

NFP 

Loop  Index  =  NFUEL  +  1 

OUTPL 

Cell  output  in  joules  in  fuel 

zones 

OUTPLOT 

Cell  output  in  joules  in  pusher  zones 

RPLOT 

Zone  radii  xn  millimeters 

TPLOT 

Total  time  in  nanoseconds 

VPLOT 

Cell  wall  velocity  in  millimeters  per 

nanosecond 
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APPENDIX  C 


Derivation  of  the  Coulomb  Cross  Sections 
for  Short  and  Long  Range  Collisions 


Figure  16  illustrates  the  hyperbolic  paths  of  two  passing 
particles  of  equal  charge  and  mass.  It  also  defines  the 
geometry  and  parameters  necessary  to  this  discussion. 
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The  asymptotes  are  in  dashed  lines,  the  impact  parameter  b 
is  defined  as  the  distance  of  closest  approach  in  the  absence 
of  other  forces,  the  angular  deflection  0  is  defined  as  the 
angle  of  deviation  from  the  normal  path,  and  the  center  of 
mass  is  point  A  .  The  angle  ^  is  an  angle  of  convenience 
between  the  asymptote  and  the  line  through  the  actual  distance  of 
closest  approach.  Note 


B  -  'TT-ZV 


(C-l) 


Extending  to  the  case  of  nonidentical  particles, 


7  tasv\  V  - 


M  b  is2. 


(02) 


where 

M  is  the  reduced  mass 

v  is  the  relative  velcoity  between  the  particles 
and  Z^  and  Z ^  are  the  charge  number  of  the  two  respective 
particles  (Ref  26:122). 

Of  interest  is  the  case  where  0  =  tt / 2  .  From  Eq  (C-l) 
it  is  apparent  that  0  =  tt/ 2  when  ^  =  tt/4  .  This  means 
tan  \p  =  1  and  therefore  the  impact  parameter  for  a  90° 
deflection,  denoted  b0  ,  is 

b  -  *  ill  ^ 

°  M  Tr*  (C-3) 


t. 
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If  one  particle  is  very  much  larger,  the  heavier  particle  may 
be  taken  as  stationary  and  the  relative  velocity  becomes  the 
velocity  of  incidence  of  the  lighter  particle.  Denoting  this 
velocity  Vj  and  noting  the  reduced  mass  under  the  condition 
m0  >>  m^  can  be  approximated  by  because 


~L  +  J-  c?  r*)t  w hen  w».» 

M, 


(C-4) 


In  other  words,  M  =  m^  .  Thus,  Eq  (C-3)  becomes 

b„= 

The  cross  section  for  such  a  collision  is  then 


(C-5) 


<77  cs 

tr  fc. 

(C-6) 

- 

v  zf  Z,V 

(C-7) 

4 

(8) 

where  E^  is  the  energy  of  the  less  massive  particle  in  the 
rest  frame  of  the  more  massive  particle. 

To  arrive  at  a  cross  section  for  a  long  range  collision, 
consider  Figure  17  illustrating  a  Coulomb  collision  in  which 
momentum  transfer  is  small  so  that  particle  1  is  basically 
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Figure  17 .  Coulomb  Collision  Assuming  Small 

Momentum  Transfer  Between  Interacting 
Particles 


undeflected  from  straight  line  motion.  This  is  essentially 
the  picture  of  an  individual  long  range  collision  between 
fast  particles  in  the  rest  frame  of  particle  2. 

Remember  the  long  range  encounter  equivalent  to  a  collective 
scatter  through  an  angle  of  90°  can  be  described  as 


( 


Ap)  =  P'NITIAI. 


(9) 


The  change  in  momentum  for  a  single  interaction,  though,  is 
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=  f  FJ* 

-4 


(C-8) 


Since  particle  1  would  apply  an  equal  force  in  both  the  x 
and  -x  directions  when  integrated  over  time,  Eq  (C-8)  can 


be  reduced  to 


•I 


F,  it 


(C-9) 


The  force  is  always  along  the  ray  r  so  that  Eq  (C-9)  can  be 


written 


- 1  Fr  A*sr>.  &  d  £ 


(C-10) 


From  Coulomb's  law,  this  can  be  expressed 


Aft  5'  Aw  a  JtL 


(C-ll) 


Note  from  Figure  17  ,  r  =  Q  ,  so  this  becomes 


j>'  X 

If  v  =  dx/dt  is  the  velocity  of  particle  1,  and  noting  that 


(C -12) 


p  =  -x  tan  0  ,  then 


b  CAC*  & 
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(C-13) 


Substituting  this  into  Eq  (C-12)  results  in 


td  = 


'V  h 


a  C*  . 

-  \  -^**1 


(C-14) 


*Ps 


2  ?,  Hz  6' 

art 


(C-15) 


In  a  large  collection  of  particles  where  N  collisions 
occur,  however,  the  total  momentum  charge  will  be 


=  Z  (Ap3)i 


(C-16) 


If  the  distribution  of  particles  is  isotropic,  &py  will 

always  be  zero  in  a  large  collection  of  particles.  The 
_  2 

quantity  (Ap^.)  ^  ,  however,  will  almost  certainly  be  non¬ 

zero  and 


(APaY  -  2 


(C-17) 


as  all  the  cross  terms  will  sum  to  zero.  This  can  be  expressed 


(apJ  =  N(*h)i 


(C-18) 


for  N  collisions.  Realizing  that  particles  can  be  incident 
from  any  angle  and  not  simply  along  the  x  axis,  this  can 
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be  extended  to 


(apY  -  N(*?)l 


(C- 19) 


so  that 

J  (ap Y  =  (Zp\  J  hi 


(C- 20) 


2  2 

and  (Ap)^  may  be  taken  as  (Ap^)  as  only  magnitudes  are 
being  considered. 

The  number  of  encounters  with  an  impact  parameter 
between  b  and  b  +  dp  per  second  depend  on  the  number 
density  of  the  interacting  test  particles  and  their  average 
energy  or  velocity.  This  will  be  described  by  a  cylindrical 
shell  of  length  A  where  X  is  the  distance  traveled  from 
one  individual  collision  to  the  next  individual  collision  or 
the  mean  free  path  or 

JN  =  2irn  <1  L  Jb  , 


Using  this  in  Eq  (C-20)  results  in 

clfZp)1  =  (£p)i  2-rrn  XhJb 

J  (5p (-  — *  bcJL 


(C- 22) 

(C- 23) 
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or 


Sir  -ZT  <el  h  a 

nj-1 


d  y 

io 


Therefore,  integration  yields 


TJ'' 


njz 


3.  yivt-A- 


xhere  A  -  XD/bmin  . 

Equation  ( C-26 )  describes  the  long  range  encounter 
equivalent  to  a  90°  scatter  when 

(zip')  - 


The  mean  free  path  is  related  to  cross  section  by 

^  '  n"<£ 


or 


(C-24) 


(C-25) 


(C-26) 


(C-27) 


(C-28) 


(C-29) 
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APPENDIX  D 


[ 


Derivation  of  the  Enhancement  Factor  for  Bosons 

A  Bose  particle  is  distinguished  from  a  Fermi  particle 
by  the  symmetry  of  its  eigenfunction.  A  Fermi  particle 
eigenfunction  is  asymmetric  and  a  Bose  particle  eigen¬ 
function  is  symmetric  (Ref  7:411).  As  a  result,  if  a  Boson 
is  introduced  to  a  Bose  distribution  of  similar  particles, 
there  is  an  enhanced  probability  the  new  Boson  will  be  born 
into  a  quantum  state  already  populated.  In  fact,  this 
increased  probability  is  highly  dependent  on  the  existing 
quantum  state  population  so  that  the  higher  the  state  popu¬ 
lation,  the  higher  the  increased  probability  or  enhancement 
factor . 

The  Bose  distribution  is  given  by  Eisberg  and  Resnick 
(Ref  7:432)  as 


I 


where  n(E^)  is  the  state  population 

E^  is  the  energy  corresponding  to  this  state 
and  T  is  the  species  temperature. 

In  a  dynamic  system  n(E^)  can  be  replaced  by  n(E^)  where 
n(E^)  is  the  average  quantum  state  population  or 
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n(£i)  =  ~ 


e"Ar-  / 


(D-2) 


Now 


-  f‘/*T 

MO  - 

I  ~  e 


CD-3) 


Thus 


n  (El)  -  n  (fr)  6 


-  ft/kT 


~£i/k  T 

e 


CD-4) 


CD-5) 


so  at  a  given  quantum  state  1 


H,  (£ii)  E,‘Ar  i 
- i 

|  +  n, (fa) 

At  another  given  quantum  state  2 

Hi  (Eii)  f“/*r  j 

- — —  e  =  1 

I  +  (£*;) 


CD-6) 


CD-7) 


Therefore 
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CD-8) 


n/e.i'l  e'V4r  r\i(En~)  e,./„ 

_ —  £  —  _ — - -  0 

I  +  n  (e,l)  I  +  ?h  (fZi) 


so  that 


therefore 


I  +  D,  (fit  ) 

nj  ( £,i) 

I  +  rT, 


CD-9) 


(f«  -£.iKt 

e 


nt  (£*;)  (  I  +  n*  (eu)) 
n,  ( i  +  nt(^)) 


(D-10) 


or 


-(Frf-ftiVfcr 

e 


7T,  (c.t)(  1  +  n>  (g*>) 
n*  (£*;)(/  *  n,  (b,l)) 


CD-11) 


Ignoring  degeneracies,  this  is  the  Boltzmann  relation 
(Ref  22:63).  Recall  the  Boltzmann  relation  gives  the  ratio 

N 

of  state  populations  for  quantum  states  characterized  by 
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w 


different  energies  and  .  Equation  (D-ll)  shows 

that  each  state  is  enhanced  by  a  factor  of  1  +  iTj  (E^)  where 
j  is  the  energy  state  to  which  it  is  being  compared.  Given 
a  situation  near  equilibrium  then,  the  enhancement  factor  is 
1  +  nj(Ej)  .  Referring  to  Eq  (D-2),  this  can  be  expressed  as 


/■+  n  (Ei)  = 


e‘/kT- 1 


+  l 


CD-12) 


or 


|  +  n  (*l)  = 


I 


ei/kT-i 


+ 


£*'At 

_e _ 2_i 

e£‘Ar- 1 


CD-13) 


which  is 


tC/hT 


1+  h  (£t)  = 


etVkr  -  i 


CD-14) 


or 

1 

j  +  h  (Ei)  ~  - 

l  ~  e 


(D-15) 


Photons  are  Bose  particles.  The  enhancement  factor  for 
photons  born  through  bremsstrahlung  or  emitted  during  a 
Compton  process  can  be  written 


|+  (Ev}  -  -e*/«r 

I  "  e 


(26) 
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where 


E 

and  T 


is  the  photon  energy  in  the  quantum  state 
is  the  radiation  temperature. 
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APPENDIX  E 


Derivation  of  Bremsstrahlung 
and  Inverse  Bremsstrahlung  Cross  Sections 

The  bremsstrahlung  cross  section  for  a  nonrelativistic 
Coulomb  collision  is  given  by  Jackson  (Ref  15:513)  as 


where 

Zj  is  the  charge  of  the  radiating  particle 

is  the  charge  of  the  particle  providing  the 
accelerating  force 

M  is  the  mass  of  the  radiating  particle 
v  is  the  velocity  of  the  radiating  particle. 

Assuming  the  radiating  particle  is  an  electron  so  that 
E^  =  Eg  ,  Zj  =  1  ,  and  M  =  mg  .  Denoting  the  charge  of  the 
particle  providing  the  accelerating  force  as  Z  and  the 
photon  energy  as  hw  =  E^  ,  Eq  (20)  becomes 

,S,)cr  Ik  ££  (_£.Yf-T—  ** 

3  Wc'/yy/fr  \  ^  )  (E-i) 

Since  the  fine  structure  constant  is  a  =  e2/ftc  ,  and  the 
classical  electron  radius  rQ  a  e2/m gc2  ,  Eq  (E-l)  can  be 
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written 


=  (E*2) 

N [  tu 

Multiplying  both  numerator  and  denominator  by  mg  and  recog¬ 
nizing 

J2_ _ \_ 

He  (E-3) 


Eq  (E-2)  becomes 

ccfcj£<)^  Ii2^  c*  ^  4  Vem£i  ")  (31) 

Ec  Ely  \  / 

if  treated  as  a  strict  equality  which  is  the  form  of  the 
bremsstrahlung  cross  section  used  in  the  text. 

The  unenhanced  bremsstrahlung  reaction  rate  is  then 

■RR^  =  r>;r>*  (E')'VV-(E',S>)  (22) 

where 

n^  is  ion  number  density 

n  is  electron  number  density 
e 

and  v  is  their  relative  velocity  or  Eg  in  the 
frame  of  the  ion. 

The  inverse  bremsstrahlung  cross  section  may  be  derived 
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using  this  reaction  rate  and  assuming  local  thermodynamic 
equilibrium.  With  this  assumption,  the  principle  of  detailed 
balance  requires  the  reaction  rates  of  bremsstrahlung  and 
inverse  bremsstrahlung  be  equal  or 

RRg  =  RRIg  (E-4) 

The  inverse  bremsstrahlung  reaction  rate  is 

K  ~  (e*) £*)  c-  0^$ (e-5) 

where 

n^(Ev)  is  the  photon  number  density 
and  cjjg  (Ee*Ev,Ev)  is  the  inverse  bremsstrahlung 

reaction  rate. 

The  ion  number  density  is  required  as  ions  are  required  to 
conserve  momentum  in  the  inverse  bremsstrahlung  process. 

The  relative  velocity  between  a  photon  and  an  electron  must 
be  c  .  Also,  as  inverse  bremsstrahlung  adds  E^  to  attain 
an  electron  at  energy  Eg  ,  the  electron  number  density 
of  interest  is  n  (E  -E  )  .  Equating  RR0  and  RRT-  , 

rtt  (fC|  E*)sr>inv(e*)h*(et~E»)c. (e-6) 


and  including  the  enhancement  factor 


n<  (S')v(E'  )  (7>  (?',£*) 

<zf(e.-E,,E,)  -  CE-7) 


For  a  higher  given  electron  energy,  this  can  be  written 

n«(EOnv(^T(V--^ 

The  electron  energies  are  in  a  Maxwellian  distribution  so  that 
the  radio  of  number  densities  can  be  expressed 

/  2  \  Ne  /  »  _  \‘A  -(E*  +  E*VkTc 


/  2  \  ^  ' 
(/5r)(fcrefi^  *  "  e 


ik.  py»  -*/bv 


nt  (£«■*£<)  _  /'£«_*_£  A*4  -Vk* 

ne(£v)  \  £<  / 


(E-9). 


(E- 10) 


At  the  new  given  energy 


<£  E'. E' )  =  t2***  *•  (^L— ) 


Substituting  these  into  Eq  (E-8) 


CE-11) 


^  <«•*>-  3  7V  r*  ) 
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The  bremsstrahlung  and  inverse  bremsstrahlung  add  or 
substract  to  the  electron  number  density  at  a  given  value 
of  electron  energy,  Eg  ,  in  four  ways.  The  number  density 
can  be  increased  by  bremsstrahlung  as  an  electron  originally 
at  a  higher  electron  energy,  Eg+E^  »  loses  energy  E^  to 
the  radiation  field.  The  number  density  can  be  increased 
by  inverse  bremsstrahlung  as  an  electron  originally  at  a 
lower  electron  energy,  Ee-E^  »  Is  boosted  to  energy  Eg 
by  a  photon  yielding  energy  E^  . 

The  electron  population  at  energy  Eg  can  be  decreased 
through  both  processes  also.  Bremsstrahlung  occurring  with 
an  electron  originally  at  energy  Eg  results  in  an  electron 
now  at  energy  Ee~Ev  •  Inverse  bremsstrahlung  occurring  with 
an  electron  originally  at  energy  Eg  results  in  an  electron 
now  at  Ee+Ev  •  The  two  processes  are  shown  schematically 
adding  and  subtracting  to  the  electron  number  density  in  a 
Maxwellian  distribution  at  a  given  energy  Eg  in  Figure  18. 

In  Figure  18  ,  processes  I  and  II  are  inverse  bremsstrahlung 
and  processes  III  and  IV  are  bremsstrahlung.  In  considering 
the  additions  or  subtractions  in  the  number  density,  the 
bremsstrahlung  and  inverse  bremsstrahlung  processes  need  only 
to  be  counted  once,  however,  if  the  distribution  function  is 
integrated  over  the  entire  range  of  Efi  . 
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Eft 


Figure  18  . 


Bremsstrahlung  and  Inverse  Bremsstrahlung 
Adding  and  Subtracting  to  the  Electron 
Number  Density 
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APPENDIX  F 


The  only  bremsstrahlung  process  considered  is  an  electron 
reacting  to  the  Coulomb  force  of  an  ion.  It  also  assumes 
the  energy  of  the  reaction  is  the  energy  of  the  electron 
relative  to  the  ion. 

To  arrive  at  the  reaction  rate  density  for  all  energies, 
the  electron  number  density  and  the  bremsstrahlung  cross  section 
must  be  integrated  over  the  entire  range  of  energies.  The 
electrons'  energies  are  described  by  a  Maxwellian  distribution. 
The  energy  range  of  particles  in  a  Maxwellian  distribution  is 
from  0  to  infinity,  which  will  be  the  limits  of  integration 
for  the  parameter  Eg  .  The  minimum  of  the  resulting  photon 
energy  is  0  for  corresponding  to  no  energy  being  yielded  to 
the  radiation  field.  The  maximum  energy  yielded  to  the 
radiation  field  is  the  entire  energy  of  the  electron  or  Ee  . 

The  limits  of  the  photon  energy  Ev  then,  are  0  to  E  , 

C 

which  will  be  the  limits  of  integration  for  the  parameter 
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The  electron  velocity  must  be  considered  in  the 
integration  of  the  parameter  Eg  .  This  can  be  written 


(F-l) 


Also,  though,  a  photon  is  a  Bose  particle,  so  if  n 
photons  already  exist  in  the  radiation  field,  the  bremsstrahlung 
process  is  enhanced  by  a  factor  of 


(26) 


Appendix  D  briefly  investigates  this  enhancement  factor. 
With  these  considerations,  Eq  (22)  becomes 


RR/5-  At  ^  ^  j 


e£^) 


(24) 


The  dimensions  are  inverse  (volume-time). 

The  total  energy  yielded  to  the  radiation  field  due  to 
bremmstrahlung  per  unit  volume  per  unit  time  can  be  attained 
by  summing  the  total  energy  of  each  photon  created.  This  can 
be  done  by  multiplying  the  reaction  rate  density  by  photon 
energy  and  integrating  over  the  range  of  possible  photon  ener¬ 
gies  or 


fa- ^  £77T7 


(> 


dE* 

S'/kTr'j 


(25) 
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Putting  the  explicit  expressions  for  the  Maxwellian  electron 
distribution  and  the  bremsstr ahlung  cross  section  into 
Eq  (25)  results  in 


Multiplying  and  dividing  twice  by  kTe  ,  substituting  for  a 
and  rQ  ,  and  reducing  yields 


Useful  simplifications  can  be  made  by  defining 

X- 


...  (29) 

(30) 

(31) 
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* 


Note  that 


W- 

i 

1/  ‘ 


E* 

ItTe 

Ev 

kTr 


(F-4) 


(F-5) 


By  substituting  these  and  further  reducing,  Eq  (F-3)  becomes 


-4  0 


d  Ev 

kJc  (F-6) 


In  preparation  to  change  variables  to  x  and  £  ,  note 

J«  -  f J  (ft) 


■'/’ J  (£) 

The  new  expression  is 

21 /3J\%'*z'(kT,y 

‘  3  \TTMj 

*  j*  0 ex&  (jr+tt) 


fa 


CF-7) 


(F-8) 


<!/  d)( 


(F-9) 
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The  upper  limit  of  the  inner  integral  is  for  the  parameter 
C  .  The  upper  limit  is  established  from  the  nature  of  the 
Maxwellian  distribution.  Remember  £  =  Eg/kTe  so  that  a 
particle  may  have  a  very  high  energy  for  a  given  temperature 
so  the  upper  limit  is  infinity.  The  lower  limit  of  the  variable 
x  is  illustrated  in  Figures  19  and  20,  graphs  of  the 
x  =  Eg/E^  versus  the  bremsstrahlung  cross  section,  so  that 
as  o0  0  -*■  x  1  (Ref  14:512). 

p 

Reordering  terms  and  substituting  results  in  the  form 
of  the  bremsstrahlung  power  density  for  a  zero  radiation 
temperature  results  in 

?«  -  n.  n.  Si  /2_f O'  m  .  c*  (kT„)* 

p  3  XW'irJ  v 


The  inverse  bremsstrahlung  reaction  rate  density  is 
stated  in  Chapter  II  as 

Jco 

n* 

o 

The  energy  limits  on  the  second  integral  differ  as  the 
photons  are  free  and  are  therefore  in  a  Planckian  distri¬ 
bution.  The  theoretical  maximum  energy  is  then  infinity. 


'CO 


nl/(£y)c^(Ee,E,)dEv 


(27) 
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Similar  to  the  bremsstrahlung  power  density,  the  inverse 


bremsstrahlung  power  density  can  be  written 


/~cz>  y»eo 

fl;  C  j  (E*)d  \  Oyf  (£<r,£v)£\>  d  S.\t 


(28) 


Substituting  explicit  expressions  into  this  equation,  it 
becomes 


This  becomes 

Tie  =  n;  n*  32.  /L2_\4zV  r.l|vie  cl  ( kT, J’’ 
r  3  Vrrmt/ 


2 

after  multiplying  numerator  and  denominator  by  E  (kTe) 

v 

Using  x  ,  £  ,  and  y  from  Eqs  (29),  (30),  and  (31),  and 

noting 


159 


(F-12) 


4-8-)= -> 


tF-13) 


*  e*V 


V 


I  -  e  w. 


A.  <7?rr  +/7)  if 

V  /-Xl  CF-14) 


The  limits  of  the  integration  over  x  can  be  seen  as 

-&!  ~  O  Y.  -  CO 


(F- 15) 


F~\j 

—  =  OD  =^>  X  -  o 


(F-16) 


Multiplying  the  second  integral  by  -1  and  reducing  yields 

1 =  n;  n<  22  (JL\-zl~<  r.1  cz  (kT*)'h 

3  \hi ^Trl 


*  Ur* 


(Vx+T  +v/x")d 


(F  -17 
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A  transformation  on  the  second  integral  may  be  performed 
by  letting 


X4  I 


then 


The  lower  limit  becomes 

X  ~  o  -*■  x'  =  I 


Now 


Realizing  x'  is  simply  a  dummy  variable  of  integration, 
Eq  (F-21)  can  be  written 


(F -18) 


(F-19) 


(F- 20) 


(F-21) 
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*  j  e‘Kx' 'V^  (/7+ y^)d; 


This  is  the  form  of  inverse  bremsstrahlung  power  density 
for  a  zero  electron  temperature.  Note  that  Eqs  (32)  and 
(33)  are  integrated  over  all  energies.  The  qualification 
concerning  the  addition  or  subtraction  of  electron  numbers 
at  a  given  electron  energy  is  therefore  satisfied. 

Equations  (32)  and  (33)  can  be  expressed  using  a  Gaunt 
factor  G(y)j  where  j  denotes  either  pure  bremsstrahlung 
or  pure  inverse  bremsstrahlung.  First,  Eq  (32)  may  be 
multiplied  and  divided  by  (kTe)*V(kTe  -  kTr)  yielding 


(F -  22) 


Note  that 


kTn-kTr 


\-x 


(F- 23) 


Therefore,  the  power  density  for  bremsstr ahlung  only 
adding  energy  to  the  radiation  field  can  be  stated 

r-‘"> c'  s  ("■-"'> 


where 


&(tV  (7 a 

Similarly,  the  power  density  for  inverse  bremsstrahlung 
only  taking  energy  from  the  radiation  field 

%  -  n; n'  ^  gsf  2V 
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where 


G  (r h  =  (v^)<U(37) 

For  practical  applications,  the  power  density  equations 
must  be  extended  from  the  limiting  cases  of  either  pure 
bremsstrahlung  or  pure  inverse  bremsstrahlung  to  a  combina¬ 
tion  of  the  competing  processes.  This  may  be  done  by  simple 
subtraction.  The  net  rate  of  energy  flow  into  the  radiation 
field  is  the  energy  that  goes  into  this  field  minus  the  energy 
that  leaves  the  field  per  unit  time.  This  is 

?,-PxP  ■  - kTr)  6  ( 4 

" Pi -  kTr)  6  <■  ^ 

(F-25) 

Denoting  the  common  coefficient  by 


C38) 
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Using  explicit  relationships  for  the  fine  structure  constant, 
the  classical  electron  radius  and  number  densities,  namely 


lie 

,0  '  Mec* 


(F- 26) 


(F -  27) 


(F- 28) 


’Em 


(F-29) 


the  bremsstrahlung  coupling  coefficient  can  be  expressed 

.  sz/z  f  e4N*, /LV!  r  G^) 

With  this  bremsstrahlung  coupling  coefficient,  the  power 
density  due  to  bremsstrahlung  processes  is 

?/"  =  (hT,- kTr) (QU)rQ (K)If) 


G(t)  =  UA-  6(4, 
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Then  referring  to  Eqs  (35)  and  (37),  this  is 


(F -31) 


Therefore 


(44) 


This  may  be  written 


6 00  -  7  l 

J*  (!-/)(,-  e  w) 


(45) 


where 


(46) 


that 

=  Ae.t  (  k-Te  -kTt)  &U) 


(43) 


APPENDIX  G 


Evaluation  of  Limiting  Values 


of  the  Bremsstrahlung  Gaunt  Factor  with  Model  Discussion 


The  Gaunt  factor  for  bremsstrahlung  may  be  stated 


u,)'  ^  H !p‘w- J  * 

At  the  limiting  value  of  y  =  0  ,  this  becomes 

e'^,(77+yrT)Jx 


which  can  be  written 


ow 

_  a 

1 <  •  fc<«  ■  * 


(Ref  11:310).  Therefore, 


6(o  )=  riA(^+^)d 


(G-l) 


(G-2) 


(G-3) 


(G-4) 


Integrating  by  parts  yields 

<00  | 


6(“'  •  Tfcr 


Letting  y  *  x  ,  this  becomes 


Another  substitution  using  Z  =  —  results  in 


■  I 


\A-z: 


d  2 


Finally,  letting  sin  0  =  x  yields 


•n/t 


G(c)-- 


A/yv\  &  C&4  e 
1  —  A****'® 


.  (*h 

G  C°)  -  \  &  de 

Jo 


so 


G(o)  -  1 
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(G-5) 


(G~6) 


(G-7) 


(G-8) 


(G-9) 


(47) 


At  the  limiting  value  of  y  =  00  ,  the  bremsstrahlung  Gaunt 


factor  is 


cO  r  00 


\{M 


(G -10) 


y°  /;  r+-  -€Uf) 


(G-ll) 


Ignoring  all  but  the  first  two  terms  of  the  series  expansion 
for  e  ,  this  can  be  written 

E-e'P-WTj 

x  J>  y*“  (G-12) 


Cancelling  the  factor  5  and  applying  L'Hospital's  rule  to 
evaluate  the  limit,  this  becomes 


g  (ce)  =  -  {V A  (v?  -v  jrr)^  ( 1  -  e* )  <1  x 

The  integration  of  g  can  then  be  done  so 

sc*)--  ^ 

/  X  ( K  -  I  \ 


(6-13) 


(G-14) 
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The  bremsstrahlung  Gaunt  factor  varies,  then,  between  the 
values  of  1  and  tt 2 / 4  .  A  numerical  integration  of  G(y) 
was  performed  using  a  two  dimensional  trapezoidal  scheme 
described  in  Reference  13.  This  data  was  fit  with  the  function 
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(126) 


Figure  21  is  a  graph  of  this  function. 

The  function  G(y)  presented  in  Eq  (44)  was  integrated 
using  a  two  dimensional  trapezoidal  integration  scheme. 

Figure  22  is  a  graph  of  the  results  of  this  integration. 

The  integration  was  designed  to  attain  a  minimal  grid 
size  at  the  areas  of  larget  G(y)  values.  As  a  result, 
areas  where  G(y)  values  were  very  small  were  neglected. 

Since  the  integral  upper  limits  are  both  ®  ,  it  is  felt  the 
neglect  of  the  aforementioned  areas  account  for  the  low  values 
of  G(y)  computed.  The  shape  is  felt  to  be  revealing, 
however,  and  it  is  this  shape  that  was  duplicated  using  the 
analytic  limits  previously  discussed. 


COMPUTED  CURVE  FOR  G(GRMMR) 


APPENDIX  H 


Derivation  of  the  Average  Energy 
Exchange  During  a  Compton  Collision 

The  average  energy  loss  of  a  photon  Eq  in  a  collision 
with  an  electron  may  be  more  easily  approached  by  first  noting 


Figure  23 ,  Symbolic  Illustration  of  Compton 
Scatter  with  Parameters 


Since  energy  must  be  conserved 

Ev  =  +  -£ 

7,m  c  (H-i) 
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and  since  momentum  must  be  conserved 


c 

and 

Eii  ~  EiL  &  +  pc+*  d> 
c,  c  ' 

Multiplying  Eqs  (H-2)  and  (H-3)  by  c  results  in 

Ev  ©  =  ^c^vi  <p 

and 

Ey  =  C+4&  +  pcc*4<f> 

Squaring  these  equations  and  then  adding  yields 

=  £*  ~2EvEi/<^©  +  Ey 

Equating  the  two  expressions  for  total  energy  gives 

Ee*r»<c*  =  (r»?c4  4  pJ-c1)/l 

By  squaring  the  two  sides  and  solving  for  p2  c2  ,  it 
be  written 


(H-2) 


(H-3) 


(H-4) 


(H-5) 


(H-6) 


(H-7) 

can 
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MICROCOPY  RESOLUTION  T  £  ST  CHART 


p*Cl  = 


(H-8) 


Ee1  +  2  Ee  roe  C1 


Since  E  =  E  -E"  this  can  be  stated 
e  v  v 


p'C2  =  ir>Z-2£y£v  4  ZM{Cl(£i;*fi)) 


(H-9) 


Equating  this  to  Eq  (H-6) , 


2  r*l  e  C*  (B»  -  £v)  r  2  Ev  Ey  0  -  -<x<  ©) 


(H-10) 


Now  6E  = 


Ev-Ev  so  that 


J' t  =  ( /  -  c&d  ©) 

We  Cz 


(H-ll) 


This  can  be  expressed 


h  (e*r.  {£}(,- 4*6) 

Me  C2 


(H-12) 


This  can  be  approximated  if  6E  is  small  by 


l/'F  -  ~~^~i  (  I  ~  lC+<i 

We  C*  v 


(H-13) 


This  approximation  gives  values  that  are  somewhat  high  as 
illustrated  by  Table  17. 
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Table  17 


Values  of  6E  Computed 
Using  Equations  (H-12)  and  (H-13) 


Ev(kev) 

6E 

Eq  13 

from 

(kev) 

6E  from 

Eq  12  (kev) 

%  Error 

1 

8.6264 

X 

10"4 

8.6189  x 

10"4 

0.09  t 

5 

2.0862 

X 

10"2 

2.0776  x 

10'2 

0.53  % 

10 

8.0668 

X 

10-2 

8.0023  x 

io-2 

0.81  1 

20 

3.1169 

X 

10"1 

3.0000  x 

10"1 

3.90  * 

The  electron  energies  used  to  compute  Table  17  are  average 
electron  energies  corresponding  to  average  scattering  angles. 
Average  scattering  angles  are  taken  from  a  graph  in  Evans 
(Ref  8:691) . 

Typical  radiation  temperatures  in  a  thermonuclear  plasma 
can  be  expected  to  be  on  the  order  of  1  kev  (Ref  9:475). 

The  error  introduced  by  Eq  (H-13)  should  not  then  be  signi¬ 
ficant  . 

Equation  (H-13)  gives  the  energy  gain  by  an  electron  at 
a  given  angle.  The  average  energy  gain  may  be  found  by 
integrating  over  angles.  The  Thomson  cross  section  has  an 
angular  dependence  of  (1  +  cos2©)  so  this  factor  must  be 
included  in  the  integration.  This  is  stated 
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rfr 

x  )  (\-c*4e>)(\  + c&4x&)  2tr^i  ©c/* 

<«*>  ■  fe  £ — — — — 

J  ((+&&&)  2ttj4a*  e  d& 


2  W|-/C^d©)(  ©c/fr 

/cff>  =  £l  _ 

'  yv>  /•! 


c*  ,  . 

\  ( I  +  £***&)  &  O 


(H-14) 


(H-15) 


The  upper  limit  of  integration  is  the  maximum  angle  of 


scatter.  Letting  cos0  *  x 


t  c  (\-*)i***, 

j  *✓-1 


^  ( i+  xl)  dx 

^  (i-X  +  Xl-X3)d* 
^  (j*f-Xa)dx 


(H-16) 


(H-17) 


<">  ■  l 


Therefore 


</£>  = 


APPENDIX  I 


Derivation  of  the 
Compton  Coupling  Coefficient 


The  reaction  rate  density  from  Chapter  II  for  Compton 

collisions  assuming  electrons  at  rest  and  photon  energies 

such  that  E  <<  m„c2  is 
v  e 


K  Rc  -  He  Dy  (Eif)  f "/h'r 

I  -  e 


(53) 


Using  the  average  energy  lost  to  the  electrons  from  the 
photon  during  a  Compton  collision  (Appendix  H)  is 


(54) 


The  power  going  to  electrons  from  photons  per  unit  volume 
is  then 


£*(*)< rTc(j£)  ,_^Av  d£, 


(55) 


No  integration  is  done  over  electron  energy  as  the  electrons 
are  assumed  to  be  at  rest.  The  Planckian  photon  energy 
distribution  is 


-I 


Sir  £v 


\  Dir  j- y  — - -  — 

s  u3  cT  e£'/kTr- i 


and  the  Thomson  cross  section  is 


(Tf  =  |Trr.z 


With  these,  the  Compton  power  density  is 


(1-1) 


*  I  -i 

u)  -  ^  f  Sir  Ev  - 8  ^  r  a.  r  JhL  * 

»«■'  n'J  ^cT  3^ 0  i-e 


~By/h'r 

^  (1-2) 


_  r«®  _  z.  ^Atr 

Pt  =  ii  -rr1  /liV  c»w  r»  \  J^e — .  <J£„ 

3  (hcl  MtC>  *  Jo  ^e'-Ar,_  i) 

Multiplying  and  dividing  by  (Mr)5  and  noting  *  =  2*h 


d-3) 


p  ii  {rfcj*.  e"~  f 

fc-  3  if  Uc/  m«cl  °J  /  t-/»n _  i}1  V 


_ (d£v 

/)*■  U"* 


Letting  y  11  Ey/^Tt  »  then 


f‘= 


d-4) 
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The  integral 


>eo 


( 


Now 


CO 


I  -  e  y 


=  2  e 


J=® 


so 


(£  *-/ 


which  is 


.  2-  co 


i  -  e j-» 


-  I(i-)e 


The  integral  then  is 

►  CD 


^  JSO  v-/© 


Denoting  the  integral  I 

X 


=  tMfeU">Va 


(1-5) 


(1-6) 


(1-7) 


d-8) 


(1-9) 


(1-10) 
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CO 

T=I 


j«o 


e(J',h[a<^p[o^] 


(I -11) 


This  is  equivalent  to 


t  i 

■»**  Jo  •»*'  4 


(1-12) 


Now 


±  -  2? 
h  44  W 


(1-13) 


(Ref  1:807),  and 

1  co 


f  e’VJ*  -  +  ' 

Jo 


d-14) 


(Ref  11:310).  Therefore  the  integral  is 


(1-15) 


Integrals  similar  to  this  are  discussed  in  Reference  23, 
Appendix  A-ll. 

The  Compton  power  density  reduces  to 
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_ - 


_ 


r.1  (kTrf 


(58) 


3  z  ^  /  i  \3  c  ne 

45  [  f)C.J 


This  refers  again  to  electrons  at  rest  or  pure  Compton 
scatter . 

The  calculation  for  inverse  Compton  scatter 
interactions  will  not  be  included  in  this  report.  Inverse 

Compton  requires  that  the  electron  have  non-zero  kinetic 
energy  prior  to  incidence  with  the  photon.  The  calculation 
first  requires  a  transformation  to  the  rest  frame  of  the 
electron.  The  energies  of  the  incident  and  emerging  photons 
are  also  transformed  in  the  electron  rest  frame.  Rates  are 
calculated  in  this  reference  frame  and  transformation  back  to 
the  original  frame  is  required.  The  result  is  (Ref  20) 


(59) 


This  equation  represents  pure  inverse  Compton  scatter. 

The  net  flow  of  energy  from  electrons  to  radiation  can 
be  found  by  subtracting  the  power  from  radiation  to  electrons 
due  to  Compton  scatter  from  the  power  from  electrons  to 
radiation  due  to  inverse  Compton  scatter  or 
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Inserting  Eqs  (58)  and  (59)  for  Pjc  and 


Pc"er  =  51  ir’  (Xfes.  r.*(hr,y fa  -  trO 

45-  )necz  N 


The  coupling  coefficient  for  Compton  scatter  then  can  be 


written 


A;f  =  ^ O.V  Cr^  ATr)+ 

Uc/  mec*  '  ' 


By  multiplying  and  dividing  by  the  Stephan-Boltzmann 


constant,  namely 


(T=  2irf 

IS  Ksc* 


(1-16) 


the  Compton  coupling  coefficient  can  be  written 


AC  -  1^8  He 


He  aIt4 


(1-17) 


Using  the  relations 


(1-18) 


N 

x  p 


(1-19) 
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APPENDIX  J 


Derivation  of  the  Lagrangian  Energy  Equation  Form 
Used  in  the  MOXNEX  Code  with  Source  Term  and  Specific  Heat  Comments 


Reference  12  gives  the  Lagrangian  form  of  the  mass, 
momentum,  and  energy  equations  as 


Oi 

T>v- 

D  t 

Vt 


♦  ±  a _c 

e  Ac 


where 


r  i  +  <L 

P*  dr 


(J-D 

(J-2) 

(J-3) 


(J-4) 


and  , 

p  is  mass  per  unit  volume 
v  is  species  velocity 
t  is  time 
r  is  length 
P  is  pressure 

E^j  is  total  species  energy  per  unit  mass. 
Multiplying  the  momentum  equation  by  v  yields 


(J-5) 


^  iz  ^  njr1-  iZ  +  Z.  -  o 
dt  dr  f>  Jr 


Expanding  the  last  derivative  in  the  energy  equation  and 
subtracting  Eq  (J-5)  results  in 


+  'viEi 


-  iz  j.  V-  -  is 


LL  +  Zbz  =  o 

H  r  dr 


(J-6) 


The  total  energy  per  unit  mass  can  be  written 


E«  =  X  *  4V': 


CJ-7) 


where  I  is  the  random  thermal  or  heat  energy.  Note  that 


c)  Er*  _ 

dt 


Ji t 


(J-8) 


d* 


ii  +  vir 

<^<t  j* 


(J-9) 


and  that 


ns  ^  d  (X4 

dr  dr "" 


(J-10) 
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(J-ll) 


Or  =  ir  j-3L  -f  or1  V  ^ 
dr  dr  dr 


Substituting  these  terms  into  Eq  (J-6)  yields 


—  +  +  2  ij£  _  0 

<i*  ir  f  &r 


If  viscosity  is  included 


6t  d-  r 


ll±  ha 
c  <JT 


Multiplying  the  mass  conservation  equation  by 
and  subtracting  gives 


P  +  q 


i2+ 

it 

■V  . 

JX 

dr 

■(^) 

if-  +i 
it 

i£ 

it 

Noting 

that 

or  = 

Jr 

d* 

this  equati 

Lon  becomes 

i£  + 

it 

dr 

H 

il 

dr 

f*  /  d* 

if 

it 

(J-12) 


(J-13) 


C J  -14) 


(  J “15) 


C J -16) 
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!£  =  lis  U 

dt  />*  it 


since  p  =  l/v.  Noting 


and 

il,  ^il+dx 

i*  <^T  it  *  <)£ 


Eq  (J-16)  may  be  expressed 


Cl  ix  4  cjX  g) V  _  ^  ^ 

it  4v  H  ~  V  7)7  eK 


(J-17) 

(J-18) 

(J-19) 


(J- 20) 


CJ-  21) 


CJ-22) 


(J-23) 
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for  an  ideal  gas,  so  Eq  (J-23)  may  be  written 

A  source  term  S  may  be  added  so  that 


This  source  term  could  be  input  from  a  laser,  a  particle 
beam,  or  another  source  of  energy. 

Also  note  that  for  an  ideal  gas 

=  12.47  Joules/deg-mole 

(Ref  23:158).  This  can  be  expressed 

cv  =  TyPart 


C J -24) 


(101) 


(J-25) 


(J-26) 


(  J -27) 
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APPENDIX  L 


Derivation  of  Equivalent  Fermi  Temperature  Equation 
for  Deuterium-Tritium  Solid  in  the  Limit  of  Full  Degeneracy 

(Te  -  0) 


The  mean  kinetic  electron  energy  in  the  full  limit 
of  degeneracy  is 


where 

Te  is  the  kinetic  electron  temperature  in  kev 
and  C£  is  the  Fermi  energy  in  the  limit  of  full 
degeneracy  given  by  Zel'dovich  and  Raizer 
(Ref  28:220)  as 


i 

%  '  ^  *  tr\th 


r\ 


% 


h  is  Planck's  constant  in  erg-sec 
k  is  Boltzmann's  constant  in  erg  per  kev 
me  is  the  electron  mass  in  grams 
ng  is  the  elctron  per  cubic  centimeter. 

Then 


(K-l) 


(K-2) 


T* 


(.S-ttldlVt 


(K-3) 


or 


dive)  ~  l‘ 46  *  /6  Ae 


(K-4) 


Now 


ne--  f  a 


~  N, 


(K-5) 


since  for  a  mix  of  hydrogen  50%  deuterium  and  50%  tritium 
A  =  2.5  and  Z  =  1  .  Therefore 


“Efc  Wwfcv.)  "  (w  *  io“)(w  x  w'5)  /■ 


V3 


(K-6) 


which  is 


^(tWsclitrz)  9  5.6S  X/O'* p2/$ 


(102) 


which  is  the  form  used  in  subroutine  HYDRO.  This  temperature 
represents  the  minimum  temperature  attainable  for  a  given 
density. 
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APPENDIX  L 


Derivation  of  the  Adiabatic 
Update  Equation  for  Radiation  Temperature 


Adiabatic  compression  or  expansion  is  described  by 


?,V,r  =  Kvj 


(L-l) 


For  an  ideal  gas,  an  equation  of  state  can  be  written 


pv*  (/-0r 


(L-2) 


where  I  is  internal  energy.  If  y  =  4/3  for  photons,  the 
adiabatic  gas  law  can  be  stated 


?,  v,  v, 1/3  =  r,v,v,,/5 


(L-3) 


and  the  equation  of  state  is 


Pv»  ii 


(L-4) 


Combining  these 


I.v/’-TaV/1 


CL-5) 


(L-6) 


v  V/j 


r- = *  (£) 


If  V2  =  +  AV  ,  this  becomes 


T,  -  I, 


\y^v/ 


(L-7) 


so  that 


— t 

+  m  J 
\/  / 


(L-8) 


Long  division  can  then  be  used  to  evaluate  the  quotient,  and 
if  &V  <<  1  ,  the  internal  energy  can  be  expressed 


(L-9) 


Again  with  the  condition  AV  <<  1  ,  the  binomial  theorem  can 

be  applied  resulting  in 


(L-10) 


The  internal  energy  is  related  to  the  radiation  temperature 


~Xr  =  ITt/V 

r  3C 


(L-ll 
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so  Eq  (L-10)  can  be  written 


which  is  equivalent  to 

and  can  be  expressed 

T  =  Tr'"')  AV/3)j  'A 


(L- 12) 


(L- 13) 


(L-14) 


(106) 


for  use  in  an  iterative  scheme  where  n  denotes  iteration 
number . 
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APPENDIX  M 
Geometry  Models  for 

Geometry  Subroutines  of  Subprogram  ALPHA1 

The  purpose  of  subroutines  XYT3  thorugh  XYT9  and 
GETA  through  GETC  in  subprogram  ALPHA1  is  to  find  the  distance 
an  alpha  particle  travels  within  a  zone.  Knowing  this  dis¬ 
tance,  energy  deposition  modelling  for  the  zone  is  done  by 
subroutines  FINDDS  and  ELOSS.  This  Appendix  will  briefly 
discuss  and  illustrate  geometry  used  in  subprogram  ALPHA1 
(Ref  16). 

The  alpha  particles  born  in  a  zone  are  assumed  to  be 
born  at  zone  center.  The  angular  paths  they  can  take  are 
illustrated  in  Figure  24.  The  angular  paths  are  most  finely 
divided  near  tt/2  . 
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Two  coordinate  systems  are  used.  One  is  for  motion 
from  the  center  of  the  zone  and  one  is  for  motion  from  a 
zone  boundary  as  it  crosses  the  boundary. 

If  the  angle  0  is  tt/2  ,  subroutine  XYT3  is  called. 
The  parameter  DIST  is  returned  as  the  distance  traveled 
in  the  zone. 


Figure  25.  Geometry  for  Subroutine  XYT3 
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If  the  angle  0  is  greater  than  tr/2  but  not  tt  , 
subroutine  XYT4  is  called.  Here  the  alpha  particle  may 
either  intercept  the  next  inner  zone  or  the  next  outer  zone. 
The  decision  is  made  on  the  value  of  the  discriminant  of 
the  solution  of  the  equation  of  the  line  of  travel  with  the 
circle  representing  the  next  inner  zone.  Subroutines  GETA, 
GETB,  and  GETC  are  called  to  compute  this  discriminant.  If 
the  particle  path  intercepts  the  next  inner  zone  wall, 
subroutine  XYT4A  is  called. 


Figure  26.  Geometry  for  Subroutine  XYT4 
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Figure  27.  Geometry  for  Subroutine  XYT4A 

If  the  next  outer  zone  wall  is  intersected,  subroutine 
XYT4B  is  called.  In  both  cases,  DIST  is  returned  as  the 
distance  traveled  in  the  zone. 

Subroutine  XYT5  is  used  if  angle  0  is  it  .  If  the 
center  of  the  microsphere  is  intercepted,  the  direction  of 
the  alpha  particle  path  is  reversed  and  the  alpha  particle 
begins  to  travel  outward.  Again,  DIST  is  the  distance 
traveled  in  the  zone. 
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Figure  28.  Geometry  for  Subroutine  XYT4B 


Figure  29.  Geometry  for  Subroutine  XYT5 
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Subroutine  XYT6  is  called  if  angle  0  is  0  .  The 


parameter  DIST  is  seen  from  Figure  30. 


Figure  30.  Geometry  for  Subroutine  XYT6 

If  angle  0  is  greater  than  0  but  less  than  tt/2  , 
subroutine  XYT7  is  used.  The  only  possibility  for  the 
path  of  the  particle  is  to  intercept  the  wall  of  the 
next  outer  zone.  The  discriminants  of  the  solution  of  the 
equation  of  the  line  of  travel  with  the  circle  representing 
the  next  outer  zone  are  required  to  compute  DIST.  Sub¬ 
routines  GETA,  GETB,  and  GETC  provide  these. 
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Figure  31.  Geometry  for  Subroutine  XYT7 

Subroutine  XYT8  is  similar  in  function  to  subroutine 
XYT4  but  is  written  for  the  case  where  subsequent  motion  of 
the  particle  is  located  on  a  zone  boundary  instead  of  a 
cell  center.  In  the  instance  the  alpha  particle  is  travel¬ 
ing  inward,  it  may  intersect  the  next  inner  zone  or  the 
next  outer  zone.  Therefore,  subroutines  XYT4A  and  XYT4B 
are  used.  Subroutine  XYT8  is  used  only  for  angles  between 
0  and  tr  . 

If  subsequent  motion  from  a  cell  boundary  is  at  the 
angle  it  ,  subroutine  XYT9  is  used.  As  with  subroutine 
XYT5,  if  the  microsphere  center  is  reached  the  path  of  the 
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Figure  32.  Geometry  for  Subroutine  XYT$ 


Figure  33.  Geometry  for  Subroutine  XYT9 


particle  is  reversed  and  the  particle  begins  to  travel 
outward.  If  the  angle  G  is  0  ,  subroutine  XYT6  returns 
the  distance  in  the  zone. 

The  subroutines  GETA,  GETB  and  GETC  return  the 
coefficients  a  ,  b  ,  and  c  from  the  quadratic  equation 


c\xx  +  bx  +  C  =  O 


(M-l) 


These  coefficients  are 


a - 


t* 


i 


(M- 2) 


(M-3) 


and 

_  x„ 

(M-4) 


where 

Tj  is  the  radius  of  zone  r^  or  r^+^ 

xQ  is  the  original  position  of  the  alpha  particlfe. 
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APPENDIX  N 


Derivation  of  the 

Electron  Temperature  Update  Equations 
The  random  thermal  energy  possessed  by  electrons  is 

Erie  =  (Ni) 

or 

kT*V  (N- 2) 

In  a  spherical  volume,  this  may  be  expressed 

Ert'--  (*-*) 

where 

r+  is  the  outside  annulus  radius 
r_  is  the  inside  annulus  radius. 

The  energy  may  be  calculated  knowing  the  radii  of  the 
volume  element,  the  electron  number  in  that  volume,  and 
the  electron  species  temperature.  In  a  non-equilibrium 
situation,  an  updated  temperature  may  be  calculated  using 
the  conservation  of  energy  expressed  in  a  flux  conserva¬ 
tive  scheme.  In  a  spherical  volume,  this  can  be  written 
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r 


+  4lT  T1  F  -4TTf;l(%At 

CN-3) 

where 

F+  is  the  outward  energy  flux 
F_  is  the  inward  energy  flux 
At  is  the  time  increment. 

Flux  is  in  energy/area/time.  Figure  34  illustrates  the 
flux  entering  and  leaving  a  spherical  volume  element. 


Figure  34 .  Flux  Entering  and  Leaving  A 
Spherical  Volume  Element 
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Fick's  law  models  flux  across  a  gradient  using  a  diffusion 
coefficient  Kg  and  may  be  stated  in  one  dimension  as 


(N-4) 


Using  this,  Eq  (N-3)  becomes 


(r;3- r_3)|nftJ*T;n 

*  «  (>■•  K.  £  |  -  £|r>* 


This  can  be  simplified  to 


CN-5) 


(N-6) 


!  "• 1  V)  ■  ^r.  (-•*  *•  %l;  r.  K.SI,.) 

CN-7) 


t;  =  r.  +  dr 
r;3-  r3  =  3t>dr 


(N-8) 


CN-9) 
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r+n-  r_3  csZr'Ar 


(N-10) 


Substituting  this  with  further  manipulations  into  Eq  (N-7), 
the  result  is 


rr-  % 


3r.^rn 


7^ 


K, 


AT« 


Arlri 


-  r-1  K< 


—I) 

Ar  lr./ 


(N-U) 


Defining 


FLUXP  =  r4l  Ke  £2*- 

*r  r. 


and 


FLUXM*  r*K.  £5- 

r. 


this  becomes 


(N- 12) 


(N-13) 


T>  Ml 

C 


Teh+  - ~ - -  /FLWXP  -  FtUXM^At 

3r.lArntk  \  ' 


(N-14) 


The  term  nfi  is  in  spherical  coordinates 


(N-15) 


which  is 
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(N- 20) 


V'-T.N 

Spitzer  gives  the  time  of  equipartition  in  Reference 

26  as 

,  _  _ 3  we  rx;  kV‘- _ /Te  + 

"  8(2t,)^n;-Hj  z!  e*x4../l«lme  "V  CN-2i) 

where 

mg  is  the  electron  mass 

m^  is  the  average  ion  mass 

n^  is  the  ion  number  density 

Zg  is  the  charge  on  the  electron 

Z^  is  the  average  charge  on  the  ion 

e  is  the  unit  of  electric  charge 

A  •  is  the  ratio 
ei 

Tg  is  the  average  electron  temperature 
is  the  average  ion  temperature. 

Since  m  <<  m.  and  assuming  T  is  not  <<  T.  ,  the 
term  T\/in^  will  contribute  insignificantly  to  the  time 
of  equipartition  and  may  be  disregarded.  Using  a  frequency 
representation  v  =  1/t  *  Eq  CN -21)  becomes 
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(N- 22) 


8  (2T r)'/*rnj/,m  ~Z  *  'Z.*  <i*  Jlti 

Ve<5  =  3m;  (  kTe  )v* 

Noting  the  two  relations 


where 

p  is  the  density 

N  is  Avagadro’s  number 
a 

A  is  the  ion  mass  number 
and  also  that  Zg  =  1  ,  Eq  (N-22)  may  be  written 


v"»=  *A'(kT,y* 


(N- 23) 


(N-24) 


(82) 


Equation  (N-20)  can  be  written 

Tc"*'=  V  +  VM  (T^-T^At  (n- 25) 

This  is  the  equation  used  to  model  equipartition  of  energy 
between  electrons  and  ions  in  subroutine  HTFLX  of  the 
MOXNEX  code. 
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t 


Derivation  of  the  Radiation  Temperature  Update  Equations 


The  energy  in  a  radiation  field  is 

Erf  -  ±?Tr4V 

C 


CO-i) 


where 


a  is  the  Stephan- Boltzmann  constant 
c  is  the  speed  of  light 
Tr  is  the  radiation  temperature 
V  is  the  containing  volume. 


The  quantity  4a/c  is  often  denoted  as  the  constant  a  or 


Erf- 


CO-2) 


In  spherical  coordinates,  this  may  be  expressed 

Er*  =  ~TT  (r+  -  r_3)  3  Tr  CO-3) 

where 

r+  is  the  outside  annulus  radius 
r_  is  the  inside  annulus  radius. 

The  energy  may  be  calculated  knowing  the  radii  of  the  volume 
element  and  the  radiation  temperature  in  the  volume.  In  a 
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non-equilibrium  situation,  an  updated  temperature  may  be 
calculated  using  the  conservation  of  energy  in  a  flux 
conservative  scheme.  For  a  spherical  volume,  this  can 
be  written 


±'rr(rS-rj)*Tc*('4l)  =  Wr«’-r»)«T,4 


where 


4  ^TT  nz  E  Ai  -  4  V  r/  £  At 


(0-4) 


F+  is  the  outward  energy  flux  at  r+  in 
energy/ area/time 

F_  is  the  negative  of  the  outward  energy 
flux  at  r_  in  energy/area/time 
At  is  the  time  increment  =  total  timen+^  -  total  time1 


Figure  35  illustrates  the  flux  entering  and  leaving  a 
spherical  volume. 

Fick's  law  models  flux  across  a  gradient  using  a 
diffusion  coefficient  Kr  and  may  be  stated  in  one 


dimension  as 


(0-5) 


1 


Figure  3S  .  Flux  Entering  and  Leaving  A 
Spherical  Volume  Element 


Using  this  equation  (0-4)  becomes 

=  i7r(r.’-r.VATr 

+  4tr (r,1 Kr—I  -r.'Kr^l 
v  Arln  Ar  lr 


Between  each  radius  increment,  the  volume  is 


I 


I 


V  -  1  tt  (  r;3  -  r_3) 


Substituting  this  with  further  manipulations  into  Eq  (0-7), 
the  result  is 


Defining 

RFi-uxp-  r/Kr 

^r|r+ 


and 

RFtWXM  =  r_*  Kr  —  I 

AT  |  r. 


this  becomes 


Tr4'M,)rTr4M4  121  (£fLU*?- K?LUm)kt 

qVc  v 


If 


(0-8) 


(0-9) 


(0-10) 


(0-11) 


(0-12) 
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V4,=  T,%A  Tr 


(0-13) 


or 


ATr  =  Tr*4 TrK 


then 


(M'^Tr4'"'+  4Tri(“’ATr 


or 


lr  -  if 


4  Tr3<'')('Tf",-Tf) 


Substituting  this  into  Eq  (0-12), 

4Tr3<"'(Tr"""-Tr1)  =  ~7  /WuXP-RFLWXM')A!t 

r  '  ^  Yt  \ 


or 


Tr  =  Tr  + - — — furmxp- rfluxm)** 

r  Vi  [ 


Since 


(0-14) 


(0-15) 


(0-16) 


(0-17) 


(0-18) 


(0-19) 
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(0-20) 


4«Tr5  =  p 


Then  Eq  (0-18)  becomes 


*  "TrV  (rflui?  -RFlcocm)  At 


(0-21) 


This  is  the  radiation  diffusion  term  of  the  radiation 
temperature  update  equation. 

Now  since 


A ec  “  Ac r  +  A<zc 


£!>  =  (t«  -Tr) 

at  cvr 


(0-22) 


the  electron-radiation  energy  transfer  portion  of  the 
radiation  temperature  update  equation  can  be  written 

T/"  -  Tr*  4  Aer  /y  "  7^)  At 

Cyr 


(0-23) 
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